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CHAPTER  1 
INTRODUCTION 

The  use  of  fiber  reinforced  composite  materials  for  structural  compo¬ 
nents  affords  the  designer  the  flexibility  of  choosing  appropriate  layer 
orientations  to  achieve  required  directional  stiffness.  These  materials 
are  particularly  attractive  in  applications  where  high  strength  to  weight 
ratios  are  important.  However,  this  flexibility  comes  an  increased  com¬ 
plexity  in  the  analysis  of  these  structures  and  numerical  techniques  are 
required. 

TI.3  goal  of  the  static  response  portion  of  the  present  project  is  the 
nonlinear  analysis  of  thin  to  moderately  thick  multilayer  composite  plate 
structures.  These  structures  may  include  cutouts  and/or  other  free  edges; 
an  important  consideration  is  the  possibility  of  severe  stress  gradients 
near  these  free  edges. 

In  order  to  obtain  accurate  prediction  of  nonlinear  structural  response, 
accurate  finite  elements,  coupled  with  an  appropriate  and  accurate  nonlinear 
analysis  scheme  will  be  needed.  In  terms  of  the  nonlinear  analysis,  alter¬ 
nate  schemes  can  be  examined  and  compared  using  problems  of  elastic-plastic 
analysis  of  single  layer  plates;  the  results  of  such  a  study  should  then 
guide  :he  selection  of  the  nonlinear  scheme  for  multilayer  plates. 

To  achieve  this  goal,  four  pilot  studies  are  described  which  will  serve 
as  building  blocks  for  a  future  static  analysis  program.  The  studies  in¬ 
clude;  (1)  the  analysis  of  edge  effects  in  laminates  under  prescribed  uni¬ 
axial  inplane  strain,  (2)  the  development  of  a  single  layer  plate  element 
with  a  straight  traction-free  edge,  (3)  the  elastic-plastic  analysis  of  single 
layer  isotropic  plates,  and  (4)  edge  singularity  analysis.  The  rationale 
for  each  of  these  tasks  is  described  in  the  following. 


2 


The  assumed-stress  hybrid  finite  element  model  is  used  in  the  first 
three  studies.  Briefly,  this  model  involves  the  selection  of  (1)  an 
intraelement  stress  field  which  satisfies  the  homogeneous  equilibrium 
equations,  and  (2)  an  intraelement  (or  element  boundary)  displacement 
field  which  yields  the  required  interelement  boundary  displacement  con¬ 
tinuity.  Elements  based  on  this  model  are  often  found  to  yield  improved 
convergence  and  intraelement  stress  predictions  in  comparison  to  corre¬ 
sponding  assiomed-displacement  elements.  The  hybrid-stress  model  is 
ideally  suited  for  multilayer  plate  applications  since  layer  stress  fields 
can  be  defined  which  satisfy  exactly  interlaminar  traction  continuity  and 
upper/lower  surface  traction- free  conditions  exactly.  In  addition  free- 
edge  conditions  can,  inprinciple,  be  satisfied  exactly. 

The  three  studies  identified  above  are  described  in  Chapters  2  through 
4  of  this  report.  Each  chapter  is  reasonably  self-contained  and  includes 
pertinent  literature  survey,  details  of  the  formulations  and  developemnts, 
example  problems,  and  summary  remarks.  The  following  is  a  brief  summary 
of  each  study,  including  the  relevance  to  the  overall  analysis  objectives. 

In  order  to  better  understand  the  nature  of  stress  distributions  in  the 
vicinity  of  cutouts  and  other  free-edges,  it  is  first  necessary  to  better 
understand  the  free-edge-vicinity  stress  distributions  in  a  more  well- 
defined  (simple)  problem.  One  such  pilot  problem  which  is  used  by  most 
investigators  is  a  multilayer  strip  (symmetrically  stacked)  of  finite 
thickness  and  width  which  are  less  than  the  strip  length.  This  is  intended 
to  simulate  a  tension  test  specimen.  Mathematically,  the  problem  is  posed 
as  a  generalized  plane  strain  analysis  (in  the  cross-section  of  the  strip) 
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in  which  the  loading  corresponds  to  a  prescribed  uniform  inplane  strain 
(normal  to  the  plane  of  analysis).  Such  a  problem  is  useful  in  the  present 
study  since  relatively  efficient  2-D  generalized  plane  strain  analysis  can 
be  performed.  In  any  finite-element  analysis  certain  of  the  elasticity 
field  equations  and  boundary  conditions  are  only  approximately  satisfied. 

The  pu?.*pose  of  one  pilot  study,  described  in  Chapter  2,  is  to  assess  the 
effects  of  enforcement  of  (1)  traction-free  edge  or  (2)  interlayer  strain 
continuity  conditions  on  the  predicted  stress  distributions  near  the  free- 
edge.  Although  the  present  study  utilizes  pseudo  2-D  elements,  the  results 
obtained  should  provide  insight  into  those  conditions  which  should  be 
incorporated  into  more  general  3-D  multilayer  plate  elements  designed  for 
use  near  free-edge  zones  of  a  structure. 

In  the  analysis  of  more  general  multilayer  plate  structures  with  free- 
edges ,  special-purpose  multilayer  plate  elements  may  be  required  along  the 
free-edges.  One  such  element  currently  envisioned  is  a  multilayer  plate 
element  for  which  the  traction-free  conditions  are  exactly  satisfied 
along  one  edge  of  the  element.  With  the  hybrid-stress  model,  this  requires 
that  a  layer  stress  field  be  defined  which  exactly  satisfies  the  traction- 
free  conditions  on  that  edge.  Numerous  plausible  stress  fields  can  be 
defined,  and  numerical  experimentation  is  required  to  identify  the  best 
stress  field.  Extensive  insight  toward  identifying  the  best  stress  field 
can  be  obtained  by  considering  a  single  layer  pure  bending  moderately 
thick  plate  element  in  which  all  components  of  stress  (bending  contributions) 
are  included.  In  Chapter  3,  an  8-node  moderately  thick  single  layer  plate 
element  with  a  straight  traction-free  edge  is  developed.  Various  stress 
fields,  each  satisfying  the  free-edge  conditions,  are  defined  and  compared 
for  selected  example  problems.  The  best  stress  field  (element)  is  identified, 
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and  can  serve  as  a  basis  for  development  of  a  special-purpose  multilayer 
plate  element  including  both  bending  and  stretching  contributions, 

The  analysis  program  envisioned  will  include  material  nonlinear 
effects.  One  phase  of  the  present  study  is  to  evaluate  alternate  nonlinear 
analysis  schemes  in  terms  of  accuracy  and  efficiency,  Alternate  hybrid- 
stress  functionals  are  defined  in  Chapter  4  for  material  nonlinear  analysis 
using  the  initial-stress  procedure.  In  this  procedure,  material  non- 
linearities  are  accounted  for  via  an  equivalent  nodal  force  vector  represent¬ 
ing  the  difference  between  an  assumed  elastic  stress  state  and  the  actual 
stress  state.  The  alternate  functionals  (approaches)  are  examined  for 
the  elastic-plastic  analysis  of  single  layer  plates,  and  the  better  approach 
is  identified.  To  extend  the  better  functional  to  multilayer  plates,  the 
single  layer  element  can  be  replaced  by  a  multilayer  element,  and  the 
elastic-plastic  nonlinear  material  model  replaced  by  appropriate  existing 
failure  and  post-failure  nonlinear  models  for  laminated  composites. 

In  addition  to  the  approaches  used  in  Chapters  2  to  4,  one  could  also 
introduce  a  special  element  at  the  free  edge  point  to  account  for  the 
exact  singular  nature  of  the  stress  there,  To  this  end,  the  first  step  is 
to  determine  analytically  the  order  of  singularity  and  the  angular  distri¬ 
bution  of  displacements  and  stresses  near  the  free-edge  point.  Although 
the  stress  singularity  analyses  of  the  free-edge  point  have  been  done  ex¬ 
tensively  for  isotropic  materials,  there  are  few  published  works  on  corre¬ 
sponding  analyses  for  anisotropic  materials.  One  of  the  difficulties  is  that 
the  Airy  stress  function  for  isotropic  materials  is  no  longer  applicable 
to  anisotropic  materials.  The  other  difficulty  is  the  possibility  of 
multiple  eigenvalues  for  the  elasticity  constants  and/or  the  order  of 
singularities.  In  Chapter  5  we  analyse  the  form  of  stress  singularities  near 
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the  free-edge  point  of  an  anisotropic  composite  wedge.  The  changes  in  the 
forms  of  stress  singularities  due  to  multiple  eigenvalues  are  presented. 

On  the  dynamic  response  of  composites,  the  emphasis  is  placed  on  the 
prediction  of  transient  wave  propagation.  Although  the  dynamic  response 
of  a  composite  subject  to  harmonic  oscillation  has  been  studied  and  well 
understood,  there  seems  to  be  no  reliable  way  to  predict  the  transient 
response  of  a  composite.  In  this  report  we  consider  two  problems  on 
transient  wave  propagation. 

In  Chapter  6,  we  present  a  theory  of  viscoelastic  analogy  for  wave 
propagation  normal  to  the  layering  of  a  finite  periodically  layered  bi¬ 
laminate.  Each  layer  of  the  bilaminate  can  be  elastic  or  viscoelastic 
materials.  The  composite  is  modeled  as  an  homogeneous  viscoelastic  material. 
The  crux  of  the  problem  is  the  determination  of  the  relaxation  function 
of  this  "equivalent"  homogeneous  viscoelastic  material.  We  obtain  the 
relaxation  function  by  comparing  the  solutions  to  the  wave  propagation 
problem  in  the  finite  layered  composite  and  in  the  homogeneous  viscoelastic 
medium.  Wave  propagation  in  the  layered  composite  is  then  obtained  by 
solving  the  wave  propagation  in  the  homogeneous  viscoelastic  medium. 

Numerical  examples  for  an  elastic  composite  show  excellent  agreements  between 
the  solution  obtained  by  this  theory  and  the  exact  solution  by  the  ray 

theorv. 

The  asymptotic  solution  in  a  serai-infinite  bilaminate  reported  in 
the  literature  shows  that  the  stress  response  oscillates  as  time  increases 
when  the  bilaminate  is  elastic.  If  the  bilaminate  is  viscoelastic,  the 
stress  response  is  monotonic.  This  presents  a  paradox  because  an  elastic 
bilaminate  is  a  special  case  of  viscoelastic  bilaminates;  and  yet  the 
asymptotic  solution  is  oscillatory  for  the  special  case  of  elastic  bi- 
laminate,  not  for  the  general  case  of  viscoelastic  bi laminate.  This  para- 
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dox  is  solved  in  Chapter  7  by  considering  the  interaction  between  the  dis¬ 
persion  and  dissipation  of  wave  propagation  in  general  viscoelastic  bi¬ 
laminates.  If  the  distance  traveled  by  the  wave  is  not  too  large,  the 
dispersion  effects  dominate  and  the  response  is  oscillatory.  As  the  waves 

travel  farther,  the  viscous  effects  eventually  prevail  and  the  response 
becomes  monotonic.  The  distance  beyond  which  wave  propagation  becomes 
monotonic  is  also  determined  in  Chapter  7. 

Finally,  in  Chapter  8  we  present  a  brief  concluding  remark  on  this 
entire  project. 
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CHAPTER  2 

FINITE-ELEMENT  STUDY  OF  EDGE-EFFECTS  IN 
LAMINATES  UNDER  INPLANE  STRAIN 


ABSTRACT 


The  hybrid-stress  finite-element  model  is  used  for  the  analysis  of  symmetric¬ 
ally-stacked  arbitrary  angle  laminates  subjected  to  a  prescribed  uniform  inplane 
strain.  The  analysis  reduces  to  a  2-D  analysis  in  the  plane  of  the  laminate  cross- 
section.  Multilayer  2-D  hybrid-stress  elements  are  developed  which  utilize  high- 
order  through-thickness  distribution  for  displacements  and  stresses.  Three  types 
of  elements  are  developed:  (I)  Standard  elements  in  which  interlayer  displacement 
and  traction  continuity,  and  upper  /  lower  surface  traction-free  conditions  are 
exactly  satisfied,  (2)  Strain  continuity  elements  in  which  the  standard  element 
is  modified  to  exactly  satisfy  appropriate  interlayer  strain  continuity,  and 
(3)  A  traction-free  edge  element  in  which  the  standard  element  is  modified  to 
exactly  satisfy  free-surface  conditions  on  a  lateral  edge  of  the  element.  These 
elements  are  applied  to  several  example  problems  and  results  are  compared  with 
existing  results  to  assess  the  effects  of  the  various  elements /strategies  used. 
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1.  Introduction 

The  cause  of  laminate  failure  in  multilayer  composite  structures  has  been  a 
source  of  speculation  for  some  time.  It  is  generally  accepted  that  high  interlami- 
neur  stress  gradients  exist  near  the  free  edge  of  these  laminates  which  are  a  cause 
or  partial  cause  of  laminate  failure.  Therefore,  although  classical  lamination  the¬ 
ory,  which  does  not  include  the  effects  of  the  interlaminar  stresses,  will  provide 
accurate  stress  predictions  away  from  the  free  edge,  it  is  unacceptable  as  a  means 
of  solution  for  this  problem  near  the  free  edge  Llj.  Consequently,  the  problem  of 
composite  laminate  failure  has  been  approached  by  several  investigators  using  a 
wide  variety  of  numerical  techniques  [l2-7j. 

Although  numerous  stress  results  have  been  presented  for  both  cross-ply  and  an¬ 
gle-ply  laminates,  not  any  set  of  solutions  has  been  accepted  as  completely  correct, 
.'lany  discrepencies  between  solutions  still  exist.  The  effects  of  approximating  known 
stress  and  strain  conditions  which  exist  along  the  surface  of  the  laminate  and  a- 
cross  interlayer  boundaries  of  the  laminate  also  remain  a  point  of  interest  when 
various  solutions  are  compared. 

In  this  study,  a  hybrid  stress  finite  element  model  is  used  to  solve  the  problem 
of  a  composite  laminate  under  uniform  inplane  strain  C^igure  1).  The  assumed  stress 
hybrid  formulation  is  well  suited  to  this  type  of  problem  due  to  the  fact  that 
stresses  and  displacements  can  be  assumed  independently  within  each  layer.  There¬ 
fore,  continuity  of  interlaminar  stresses  and  continuity  of  inplane  strains  across 
interlayer  boundaries  as  well  as  continuity  of  displacements  need  not  be  approxi¬ 
mated,  but  can  instead  be  exactly  satisfied.  In  addition,  the  traction  free  edge 
condition  along,  the  outer  surface  of  the  laminate  can  be  exactly  enforced. 

A  special  purpose  multilayer  element  is  developed  which  satisfies  continuity  of 
interlaminar  stresses  across  interlayer  boundaries  and  which  also  satisfies  the 
traction  free  edge  condition  along  the  top  and  bottom  of  the  laminate.  A  finite 


(a)  LAMINATE  SUBJECT  TO  UNIFORM  INPLANE  STRAIN  fix 


.  TRACTION 
FREE  EDGE 
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(b)  PLANE  OF  ANALYSIS 


Fig.  1  Problem  Definition 
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element  mesh  composed  of  these  multilayer  elements  comprises  the  first  approach 
used  to  solve  the  above  stated  problem-  Two  additional  analyses  are  also  considered. 
The  first,  modifies  the  multilayer  element  to  include  the  traction  free  edge  condi¬ 
tion  along  the  lateral  side  of  the  laminate.  This  element  is  then  coupled  with  the 
standard  multilayer  element  to  comprise  the  second  approach  to  the  problem.  In  the 
third  analysis,  the  origiaal  (standard)  multilayer  element  is  modified  to  include 
the  continuity  of  inplane  strains  (computed  from  stresses)  across  interlayer  boun¬ 
daries.  A  mesh  composed  of  only  these  strain  continuity  elements  comprises  the  third 
approach . 

In  the  chapters  to  follow,  the  formulation  of  the  special  purpose  multilayer  ele¬ 
ment,  including  all  additional  modifications,  is  presented.  Tost  cases,  including 
cross -ply  and  angle-ply  examples,  are  described  and  stress  results  are  shown.  The 
differences  and  similarities  between  the  three  approaches  are  discussed.  Observa¬ 
tions  between  the  stress  results  for  different  laminate  stacking  sequences  are  also 
discussed  where  appropriate. 
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2,  Formulation 

A.  Energy  Functional  and  Element  Matrix  Definitipns 

The  general  formulation  of  the  hybrid-stress  model  for  cross-ply  laminates  sub¬ 
ject  to  uniform  inplane  strain  is  given  in  ^9,9].  That  development  will  be  expanded 
here  to  include  the  more  general  angle-ply  laminate  case.  The  enforcement  of  cer¬ 
tain  displacement  and  stress/strain  continuities  will  also  be  incorporated  into  the 
formulation. 

The  basis  of  the  assumed-stress  hybrid  formulation  is  a  modified  complementary 
energy  principle.  This  is  a  two  field  principle  for  which  intraelement  equilibra¬ 
ting  stresses  and  interelement  compatible  displacements  can  be  assumed  independent¬ 
ly  within  each  element.  In  the  case  of  multilayer  structures  which  are  subdivided 
such  that  a  single  element  is  composed  of  a  number  of  layers,  the  energy  must  be 
simimed  over  the  numbers  of  layers  within  each  element  as  well  as  the  number  of  ele¬ 
ments.  Thus,  the  hybrid-stress  functional  for  multilayer  structures,  disregarding 
external  forces,  is  given  as 


CD 


z 


(r‘ 


=  the  siommation  over  the  number  of  elements 
=  the  summation  over  the  number  of  layers 

■  the  volume  of  the  i  layer  for  the  n  element 

»  the  components  of  the  stress  vector  for  the  i^^  layer 

a  the  components  of  the  strain  vector  in  the  i^^  layer  based  on 

displacements 

til 

a  the  material  property  matrix  for  the  i  ”  layer  in  the  global 
x,y.z  system 
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The  effect  of  the  uniform  inplane  strain,  c^,  will  be  incorporated  directly 
into  the  complementary  energy  functional.  In  order  to  accomplish  this,  first  con¬ 
sider  the  expression  for  the  material  property  matrix,  S.  A  material  axis  system 
is  defined  such  that  Xj  is  in  the  fiber  direction  for  the  layer  and  X2  is  perpen¬ 
dicular  to  the  fibers  (inplane).  Then  the  stress-strain  relations  in  this  material 
.system  are  ClO]: 


In  general,  fibers  are  not  oriented  in  the  global  coordinate  direction.  The 
material  property  matrix  in  terms  of  stress/strain  in  the  global,  x,y,z,  system 
(element  coordinate  system)  can  be  defined  by  applying  the  appropriate  transforma¬ 
tion  laws,  and  is  given  by  ClOl: 


where  the  S's  are  a  function  of  the  original  terms  in  S  and  the  layer  fiber  orien¬ 
tation.  The  terms  in  S  are  defined  in 
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Imposing  the  condition  that  the  inplane  strain  calculated  from  stresses  must 
equal  the  prescribed  inplane  strain  permits  elimination  of  according  to: 


■  z  ^  ’  t  '  t: 


In  matrix  form  this  can  be  expressed  as: 


Substituting  this  relationship  into  the  first  term  in  (1)  leads  to: 


or  simplifying: 
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Using  equation  (5)  and  substituting  the  linear  strain  displacement  relationships, 
equation  (1)  becomes: 


I  I  AC 


,  1 


L  <7^  <J\  <7^  J 


+  fi,  ]  L<r7 


<r^r 


IhV 


I- 

S 

T,. 

I± 

X.  j 


<^1 

<r^, 

<r,y  J 


dA 


(8) 


-L 


l_Oy  0% 


d» 


<iA 


] 


9w  3v 

^1 

where  the  linear  strain-displacement  relations  are  based  on  a  displacement  field 
corresponding  to  uniform  inplane  strain  in  the  x  direction  Cil!: 


u"  Cx,y^i)  -  x?<  +  u‘  (y.^) 
v‘  (x,  /,  2)  =  v‘  ^y,0 
w‘  (x,y,i)  -  w‘ 


(9) 


Note  that  the  constant  terms  have  been  dropped  in  view  of  the  fact  that  it  is  the 

variation  of  II  ^  which  is  ultimately  of  interest, 
me 

The  stresses  are  now  expressed  within  each  layer  in  terms  of  a  set  of  stress 
parameters,  S,  such  that  the  homogeneous  equilibrium  equations  are  exactly  satis¬ 
fied.  The  reduced  equilibrium  equations,  where  the  stress  components  are  indepen¬ 
dent  of  the  X  direction  (in  view  of  equations  (9))i  are: 
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£Sy  +  in*  -  o 

££y  +  ,  o 

The  stresses  are  therefore  given  as: 


(10) 


(11) 


where  the  form  of  the  matrix  P  is  chosen  so  that  equilibrium  is  satisfied. 

The  displacements,  also  defined  independently  within  each  layer,  are  given  in 
terms  of  nodal  displacanents,  q: 


4 


I ‘ 


(12) 


a*  •  A/ 

The  displacement  interpolation,  N,  is  formed  such  that  continuity  between  elements 
is  guaranteed.  Using  the  strain  displacement  relations,  the  matrix  B  is  defined  as: 

il 

ay 

»• 


H 


2  i. 


(13) 


Substituting  (11)  and  (13)  into  (8)  leads  to '• 


TT. 


ilA 


(14) 
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where  represents  the  first,  second  and  fifth  rows  of  P^. 

Simplifying  this  expression  by  defining  the  following  layer  matrices^ 


tL'  - {  i‘V  r  ■ 

"lu 

(1)^ 


(15a) 


(15b) 


(ISc) 


leads  to: 


Up  to  this  point,  amd  3^  are  independent  from  layer  to  layer.  Recalling  that 
a  single  element  is  composed  of  a  number  of  layers,  improved  results  can  be  found 
by  enforcing  certain  continuity  conditions  along  interlayer  surfaces  such  as  con¬ 
tinuity  of  displacements,  continuity  of  tractions,  and  continuity  of  strains  from 
stresses.  Appropriate  relationships  between  the  stress  parameters,  3^,  and  nodal 
displacements,  q^,  for  a  layer  and  stress  parameters,  3j^,  and  nodal  displacements, 
q  for  the  laminate  can  be  defined  in  terms  of  these  continuity  conditions.  For  the 
present  case  of  prescribed  1  they  can  be  expressed  in  the  following  form: 


(17a) 
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(17b) 


Substituting  C17a)  and  (I7b)  into  (16)  and  defining  the  following  laminate  matrices: 


H  .r  • 

t  _  1. 


H  .  Z  T,  H  T, 

^ 


(18a) 


^  ^  "H  ^ 

imi  ^ 


r  .  T 

‘  v-  *  /*  ». 


T,  -  ?  7^‘ 


(18b) 


(18c) 


(18d) 


leads  to: 


TT, 


«e 


f.  /  a  &■  5 


(19) 


Recalling  that  the  3's  are  independent  from  element  to  element,  the  stationary 
condition  of  with  respect  to  eliminates  the  stress  parameters  on  an  element 
level.  Therefore,  for  an  arbitrary  53  ft  0: 


A  *  ^  P-  ^  tc‘  £* 

Substituting  this  back  into  equation  (19)  yields: 


where : 

9* 


(21b) 
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define  the  element  stiffness  matrix  and  the  element  load  vector  (due  to  e  )  re- 
spectively. 

After  assembly  of  all  elements,  the  solution  of  equation  (21)  results  in  the 
values  of  the  nodal  displacements.  Stresses  can  be  found  by  using  (20)  to  obtain 
the  laminate  3's,  (17a)  to  obtain  the  layer  S's,  and  finally  using  (11)  for  the 
stress  field. 

B.  Displacement  Interpolation  and  Enforced  Displacement  Continuities 

The  same  high  order  through  thickness  displacements  assumption  used  in  [Sj  and 
for  the  displacements  v  and  w  (in  the  y  and  z  directions  respectively)  is  util¬ 
ized  in  this  study.  Furthermore,  the  u  displacement  (in  the  x  direction)  is  chosen 
to  be  of  the  same  form  as  the  v  displacement.  The  position  of  the  nodal  displace¬ 
ments  for  each  layer  is  shown  in  figure  2b.  Note  in  figure  2a  that  the  total  height 
of  the  laminated  element  is  H,  and  the  total  length  is  1.  Also  shown  is  the  z  co¬ 
ordinate  for  the  bottom  of  layer  i  which  is  designated  and  that  for  the  top  as 
^i+l‘  convenience,  a  local  codnrditiate,  C,  is  defined  whose  origin  is  at  the  mid- 
svirface  of  each  layer.  The  local  system  is  defined  such  that; 


where 


;  =  1 


2  =  h, 
1 

z  =  h. 


The  displacement  assumption  in  terms  of  the  local  coordinate  C,  where  the  dis- 

3  2 

placements  u  and  v  are  of  order  z  and  the  displacement  at  w  is  or  order  z'",  is 

given  as; 


trUf)’  it  1T‘ 

^  (<!*  zif  -  ir  {■!  ‘7  j/  -  %  ] 


(23) 
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(b)  ELEMENT  DEGREES  OF  FREEDOM  FOR  A  TYPICAL  LAYER 


Fig.  2  Element  Geometry  and  Degrees  of  Freedom 
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*7!:  [(-I  *f  *  9f‘ 

w‘(J)  K.‘  *('(-/')k‘  •'  T  ^'-’7/3 


where  the  u  displacement,  as  previously  noted,  is  of  the  same  form  as  the  v  dis¬ 
placement. 

This  displacement  interpolation  insures  compatibility  between  elements.  To  en¬ 
force  displacement  continuity  between  layers,  however,  the  nodal  displacements  for 
layer  i  must  be  related  to  those  of  layer  i+1.  As  can  be  seen  from  figures  2a  and 
2b,  this  would  require: 


(  ,  Vi  ,  Wy)  =  f  U,  ,  V,  ,  W,  ) 


This  is  the  expression  which  relates  the  layer  degrees  of  freedom  to  the  laminate 
degrees  of  freedom  and  it  defines  the  transformation  matrix  T^  used  in  (17b) .  When 
the  nodal  displacements  are  given  the  ordering, 

*  if  )  (2S, 


where  q^  represents  the  displacements  at  the  bottom  node  of  the  layer  and  re¬ 
presents  the  displacements  at  the  upper  most  node  of  the  layer,  the  matrix  T^  be¬ 
comes  a  Boolean  matrix.  The  products  defined  in  (18b)  and  (ISd)  can  therefore  be 
accomplished  by  the  use  of  assembly  pointers  which  simply  position  the  layer  contri¬ 


butions  into  the  appropriate  laminate  positions  within  G  and  J^. 
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It  should  be  noted  that  the  displacement  interpolation  is  written  in  terms  of 
the  normalized  coordinate,  ^,and  consequently  the  matrix  B,  which  is  defined  by 
this  interpolation  using  the  strain  displacement  relations,  is  also  in  terms  of 
i;.  In  order  to  perform  the  integral  (ISb),  a  transformation  of  coordinates  must 
be  defined.  Furthermore,  since  all  of  the  integrals  will  be  evaluated  by  use  of  a 
Gauss  integration  scheme,  the  limits  of  the  integral  must  be  from  -1  to  1.  There¬ 
fore,  an  additional  local  coordinate*  s,  is  defined  in  the  y  direction  such  that 


y  -  2^2^ 

The  Jacobian  of  the  transformation  is  then  given  as: 

-X  (k  -Ki*.) 

4 


C26) 


C27) 


Thus,  the  B  and  T^  matrices  have  been  defined.  It  remains  to  define  the  stress 

assumption  and  the  continuity  requirements  placed  on  stress  and  strain. 

C.  Stress  Assumption  and  Enforced  Traction  Continuities 

The  basic  stress  assumption  used  in  CSl  and  C93  have  been  employed  here  where 

3  5  4 

Oy,  a^,  are  of  order  z  ,  z  and  z  respectively.  The  additional  stress  compo¬ 
nents  used  here,  namely  a  and  a  are  chosen  to  be  in  the  same  form  as  o  and 
’  ^  xy  xz  y 


yz 


respectively.  Although  the  basic  stress  assumption  is  the  same,  it  has  been 


recast  in  the  following  form  for  reasons  which  will  be  discussed  subsequently.  Also, 
recall  that  the  following  stress  assxanption  is  for  a  typical  layer  i: 


hi  - /,  "")(l  *  ITi  ’  J  C28a) 
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+  r h-t)  *  -  t/.‘  y 

-f/tV  in-S)U-3t)*T(k'y‘  y’)(i-f)(-J*zX*s-f') 

+  ^  y'  */4r 

-  (js  1  i)x(i--f)(i+x) 

^  /+  f’)(s’-3'f)  *  'iiCjSi'’  * jSr  y  '>  C  l-i)' ( t  -I-  ?  X*  3f'  ) 

’■*  lySi""  *^“'y)(H-Xy  it-yx+sX')*  y) 

Cl-X)'(I*X)‘*T  (/■,‘  ^3^'y)(l.f)‘-(,^X)’ 

«»'=■  T^‘  (I-X)'(Z*X)*  ifi."'li*X)'(z-f)--h.(^y*^,‘y') 
(l-f)'(l*SX)(l->3X)-kL^,''y^//y‘)(uX)‘(7-3f)(/-3fi 

*(^'  -^‘"")y^(l-X)'(l*X]'- 

-  0-X)'(l*X)  (i-X)(i+X)‘- 

*  *  X^/y')  0-f)'(n-X)(/trX)-^  (l^^'y  *3^^y^ 

Ci-i)( /*/;•  f/- r/J  -f^J{/-xr(ux)- 
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*  h  ^,‘-^‘*')yV/-i-)  *  rfiUl-f) 


(28d) 


^(/,1  /+A  V‘  zf-/) 

^  c/r'  ^^:y^)(f’-f) 


t  ^ 


(r„' n  T  (S^*/s,y  *js,'r)(f’-3f*i.') 

(28e) 

*  i  ('yS'*'  *  ^‘7  */;•'> 

*  *5/;'^  *3^  y‘)(.J^<-f‘^f-,) 

+  ^  *3y6u^  y')i-f^-f''*  J-h  1) 

*3/3,^ y')(-f*^2.r-l) 


where  for  convenience  is  defined  as  the  half  thickness  of  layer  i.  Also,  the 
stresses  defined  above  exactly  satisfy  the  equilibrium  equations  (equations  (10]). 


24 


The  special  form  into  which  the  stress  assumption  has  been  cast  easily  allows 
the  identification  of  all  the  stresses  at  the  interlayer  boundaries.  This  greatly 
simplifies  the  enforcement  of  given  stress/strain  continuities  along  the  layer  in 
terfaces.  The  interlaminar  transverse  shear  and  normal  stresses  axe  only  related 
to  appropriate  S’  terms  at  the  interlayer  boundaries. 


kr  2-Ai  ot  y*  “/ 


Hr  I 


tfy,'  ■/,  /  */*  y‘ 

•  3" 


/* 


(29) 


Similarly,  the  inplane  stresses  are  only  related  to  appropriate  S  terms  at  the  in- 


terlayer  boundaries. 

kx  a  •  A; 

4r  2  - 

0£ 

y»/ 

Oy  y  */St  y  y 

^  A  ft  ^  ^  ft 


Vh  / •'A  y y 

In  this  section,  only  the  enforcement  of  interlaminar  stress  continuities  and  sa¬ 
tisfaction  of  traction  free  conditions  along  the  top  and  bottom  of  the  laminate 
will  be  considered.  The  purpose  for  casting  a  and  a  in  the  above  form  at  the 

y  ^y 

layer  interfaces  will  be  discussed  later. 

The  traction  free  conditions  along  the  top  and  bottom  of  the  laminate  are  easi¬ 


ly  satisfied  as  a  result  of  the  form  of  the  given  stress  assumption  (28) .  Referring 
to  equations  (29)  these  conditions,  along  the  top  of  the  laminate,  result  in: 
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(Ty*'’  )  -0 

—  tUi  —  —  M*  1 

yf,  ^  o 

^  »  0 

(V  »  Avto)  *0 

3"^'  — Ato#  -/.to-/ 

V'  'A 

(31a) 


where  N  is  the  total  number  of  layers.  Enforcing  the  same  condition  along  the 
bottom  of  the  laminate  leads  to; 


/^/  • ^  “A  ’A  'A  “A  *  ^ 


(31b) 


This  is  accomplished  by  not  assembling  the  contribution  of  these  3's  into  the 
H  and  G  matrices  by  setting  the  appropriate  assembly  pointers  equal  to  zero.  Thus, 
the  traction  free  conditions  along  the  top  and  bottom  of  the  laminate  are  exactly 
satisfied. 

Continuity  of  tractions  across  interlayer  boundaries  (i.e.  a  ,  a  ,  a  )  can 

to 

be  satisfied  exactly  by  again  considering  equations  (29).  If  the  stress  parameters 
are  ordered  in  the  following  manner: 


where  the  form  of  P  is  defined  by  the  stress  assumption,  then  the  transformation 
matrix  T^  is  a  Boolean  matrix  and  the  matrix  is  zero  (see  equation  (17a)).  A- 
gain,  since  the  transformation  matrix  is  Boolean,  all  matrix  multiplications  in¬ 
volving  the  matrix  T^  can  be  accomplished,  as  previously  noted,  by  the  use  of  as 
sembly  pointers. 
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Also,  as  previously  discussed,  all  integrals  involving  the  P  matrix  will  be  e- 
valuated  using  a  Gauss  integration  scheme,  and  therefore  use  will  be  made  of  an 
additional  local  coordinate,  s,  and  the  Jacobian  of  the  transformation. 

In  general,  therefore,  the  complete  formulation  of  a  multilayer  structure  sub¬ 
ject  to  defined  displacement  and  traction  continuities  across  interlayer  botmdar- 
ies  has  been  described.  Besides  the  conditions  already  satisfied,  however,  use 
can  be  made  of  additional  stress/strain  conditions  which  can  be  imposed  on  the 
structixre  to  perhaps  provide  a  more  realistic  solution  to  the  problem.  Two  of 
these  more  specialized  conditions  are  the  enforcement  of  a  traction  free  edge 
condition  along  the  right  and  left  edges  of  the  laminate  (Pigure  1),  and  the  en¬ 
forcement  of  certain  strain  continuities  across  interlayer  boundaries.  These  ideas 
will  now  be  discussed  in  more  detail. 

D.  The  Traction  Free  Edge  Condition 

The  traction  free  edge  condition  along  the  lines  y»b  and  y»-b  as  shown  in  fig¬ 
ure  C3a) ,  is  enforced  in  a  manner  similar  to  that  of  the  traction  free  edge  con¬ 
ditions  along  the  top  and  bottom  of  the  laminate.  The  stresses  to  be  considered 
along  yab  and  y»-b,  however,  are  not  of  a  special  form  as  was  the  case  along  the 
top  and  bottom  of  the  laminate.  The  y  locations  of  the  traction  free  edges  must, 
therefore,  be  substituted  into  the  general  stress  assumption  (28)  In  order  to  de¬ 
fine  the  S's  which  must  be  set  equal  to  zero  in  order  to  enforce  the  conditions. 

To  facilitate  this  procedure,  a  local  y'  coordinate  is  adopted  at  the  left  edge 
of  the  lamiaate.  Also,  due  to  the  symmetry  of  the  laminate  about  both  the  y  and  z 
axis,  only  the  upper  left  hand  portion  of  the  laminate  need  be  analyzed.  There¬ 
fore,  as  shown  in  figure  (3b),  the  line  y»0  is  the  only  traction  free  edge  which 
need  be  considered.  As  a  consequence  of  the  translational  invariance  of  the  ele¬ 
ment  being  used,  the  stress  assumptions  are  written  in  terms  of  y  and  the  trac¬ 
tion  free  edge  condition  is  applied  to  this  set  of  equations  resulting  in: 


(a)  Pf^OBLEM  DEFINITION:  SYMMETRIC  LAMINATE 
UNDER  UNIFORM  Ex 


(b)  FINITE  ELEMENT  MODEL  OF  THE  UPPER  LEFT-HAND 
PLANE  TO  BE  ANALYZED 


Fig.  3  Finite  Element  Model 
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Thus,  if  the  elements  are  numbered  from  left  to  right  (figure  3b),  the  contribu- 
tion  of  the  above  S's  (33)  into  the  H  and  G  matrices  is  not  taken  into  account  for 
the  first  element.  Again,  this  is  done  through  the  use  of  assembly  pointers.  All 
other  features  of  the  element  remain  unchanged,  therefore  this  traction  free  edge 
(TFE)  element  remains  comparable  with  the  other  elements  to  be  used  in  the  mesh. 

One  additional  consideration  is  of  importance  when  using  the  TFE  element.  The 
work  done  by  the  tractions  is  given  by: 

as  it  appears  in  equation  (19).  For  the  case  of  the  TFE  element,  where  no  work  is 
done  along  the  left  edge,  (34)  takes  on  the  special  form: 

1  2 

where  ^  and  represent  the  displacements  of  nodes  1  and  2  respectively.  Substi¬ 
tuting  the  above  form  of  the  G  matrix  into  the  expressions  for  the  element  stiff¬ 
ness  matrix  and  the  element  load  vector  leads  to: 


(36) 


Recalling  that  the  first  element  is  taken  to  be  the  TFE  element,  the  general  form 
of  the  system  of  equations  after  assembly  is: 
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C37) 


where  K*,  q*,  and  Q*  represent  the  assembled  stiffness  matrix,  displacement  vec¬ 
tor,  and  load  vector  after  the  contributions  due  to  node  1  are  removed. 

In  practice,  this  system  of  equations  can  be  solved  in  the  conventional  manner 
by  setting  the  displacements  at  node  1  equal  to  zero  and  setting  the  diagonal  stiff¬ 
ness  terms  equal  to  unity.  This  results  in  the  system  of  non-singular  equations: 


(38) 


.\s  can  be  seen  from  equation  (20),  the  stress  parameters  will  also  be  independent 
of  the  displacements  at  node  1: 


-f£  J  f  f j  )  -  if' 5  (3. 

The  calculation  of  the  stresses  (11)  which  are  the  quantities  of  interest, 
will  not  be  affected  by  artificially  constraining  the  degrees  of  freedom  at  node 
1.  Thus,  the  singularity  which  occurs  in  the  assembled  stiffness  matrix,  as  a  re¬ 
sult  of  the  enforcement  of  the  traction  free  edge  condition,  can  be  eliminated  in 
this  way.  Solutions  for  q^  can  subsequently  be  obtained,  but  they  are  not  calcu¬ 
lated  or  used  in  this  analysis. 

E.  The  Strain  Continuity  Condition 

In  addition  to  the  stress  continuities  previously  defined  across  layer  inter¬ 
faces,  certain  strain  continuities  which  are  known  to  exist,  can  also  be  exactly 
enforced.  The  layers  of  the  laminate  are  assumed  to  be  perfectly  bonded  and  there 
fore,  the  u  and  v  displacements  along  any  xy  plane  (see  Figure  1)  are  smooth  con- 


i 


30 


tinuous  functions.  If  this  is  the  case,  their  derivatives  with  respect  to  x  and  y 

are  also  continuoiis.  Applying  this  to  the  linear  strain-displacement  relations, 

based  on  the  displacement  field  given  in  (9) ,  shows  that  both  e  and  v  are  con- 

y  xy 

tinuous  across  interlayer  boundaries. 

Using  equations  (3)  and  (4) ,  and  the  definition  of  the  matrix  R  given  in  (6)  and 

(7),  leads  to  the  following  expressions  for  e  and  y  (from  stresses) : 

y  xy 

^  - 

<ri  (TJy  (40) 

Continuity  requires  (see  Figure  2a) ; 


(2  •hi  )  -  6/’* 
/“jy  Va  -  A;  )  -  t-xy 


Substituting  (40)  into  (41)  and  recalling  that  a  is  continuous  across  interlayer 

boundaries  (i.e.  )  leads  to: 

z  z  z 


-g.'"  ay  .  zj'  <nr  ‘  ‘ 

)  - 

'  Si*  ii'"*  »  O 

■'■(x*  ~  J  *  o 


(42) 


where  all  stresses  are  evaluated  at  2  =  h. 

1 

Further,  substituting  (30)  into  (42) ,  and  collecting  coefficients  in  y,  results 
in  four  sets  of  two  simultaneous  equations  in  t's  and  J's.  Solving  each  set  of  e- 
quations,  eight  of  the  ^'s  are  expressed  in  terms  of  the  remaining  ^’s  and  T's. 
This  results  in  the  following  relationships  between  betas  of  layer  i  and  i-1 
required  to  satisfy  the  strain  continuity  condition: 


/\ 
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The  laminate  betas  correspond  to  the  set  of  stress  parameters  remaining  after  the 
interlayer  traction  and  strain  continuity  conditions  have  been  satisfied.  Equa¬ 
tions  (43)  define  the  transformation  matrices  and  shown  in  (17a).  Again, 
the  calculations  given  in  (18)  are  accomplished  both  by  use  of  the  transforma¬ 
tion  matrices  and  by  the  use  of  assembly  pointers.  All  displacement  calculations 
are  now  done  with  respect  to  a  complete  set  of  laminate  betas. 


Recall  that  stresses  however  are  assumed  independently  within  each  layer  and, 

therefore,  ntust  be  calculated  independently  for  each  layer.  This  requires  that  the 

original  set  of  layer  betas  (for  each  layer  where  stresses  are  calculated)  must  be 

calculated  using  (17a)  where  all  quantities  on  the  right  hand  side  are  known. 

Prior  to  the  consideration  of  the  strain  continuity  condition,  the  transformation 

matrix  T^  is  a  Boolean  matrix  and  all  terms  in  are  zero.  When  this  is  the  case, 
-s  ~s 

the  actual  values  of  the  layer  betas  remain  unchanged,  and  it  is  only  a  matter  of 
identifying  those  betas  in  the  laminate  set  of  betas  which  correspond  to  each  sub¬ 
sequent  layer  by  use,  of  the  assembly  pointers.  In  the  case  where  strain  continutiy 
is  enforced,  however,  ^  and  3g_j2  have  been  redefined  in  terms  of  the  set  of 
laminate  betas.  The  relationship  given  in  (17a)  must,  therefore,  be  used  in  con¬ 
junction  With  the  assembly  pointers  to  recalculate  the  layer  betas  to  be  used  in 
the  stress  calculation.  After  the  layer  betas  have  been  identified,  the  calcula¬ 
tion  of  stresses  is  identical  to  that  previously  described  (11). 


Exanrile  Problems  and  Numerical  Results 


Three  basic  approaches  can  be  defined  for  the  analysis  of  cross-ply  and  angle- 
ply  laminates.  Each  approach  exactly  satisfies  traction  free  edge  conditions  along 
the  top  and  bottom  of  the  laminate  and  continuity  of  displacement  and  interlami¬ 
nar  stresses  across  layer  interfaces.  The  element  which  satisfies  these  conditions 
will  be  termed  the  standard  element.  Other  special  stress  and  strain  conditions 
are  then  imposed  in  addition  to  the  above  continuities  and  traction  free  conditions. 
The  three  analysis  used  are:  1]  MO  SPECIAL  ELEMENT  (MSE)  analysis  where  the  stan¬ 
dard  element  is  used  throughout  the  mesh,  2)  TRACTION  FREE  EDGE  (TFE)  analysis 
where  one  special  element  exactly  satisfying  free  edge  conditions  along  the  side 
of  the  laminate  is  used  in  conjunction  with  the  standard  element  throughout  the  re¬ 
mainder  of  the  mes,h,  and  STRAIN  CONTINUITY  (SC)  analysis  where  th^  continuity 
of  inplane  strains  across  interlayer  boundaries  is  exactly  satisfied  for  all  ele¬ 
ments  . 

Five  test  cases  have  been  chosen  to  illustrate  the  effects  of  exactly  satisfy¬ 
ing  different  stress  and  strain  conditions.  Due  to  the  symmetry  of  the  plane  of 
analysis,  only  one  quarter  of  the  cross  section  need  be  modeled.  As  previously 
stated,  the  upper  left  hand  plane  is  analyzed  in  order  to  facilitate  imposing  the 
traction  free  edge  condition.  To  be  consistent  with  available  results  on  this  top¬ 
ic,  results  corresponding  to  the  upper  right  hand  portion  of  the  cross  section  are 
plotted.  Thus,  the  center  line  is  along  y=0  and  the  traction  free  edge  is  along 
the  line  y»b  (Figure  3a).  The  test  cases,  shown  in  figure  4,  are  the  4-layer  cross- 
plys  CSO/OJ^  and  CO/SO]^,  the  4-layer  angle-plys  L4S/-45]g  and  the  8-layer  lami¬ 
nates  l90/0/-4S/45]  and  II45/-4S/0/90]  . 

s  s 

The  finite  element  model  is  shown  in  figure  3b.  The  plane  of  analysis  is  broken 
into  two  regions.  More  elements  of  a  smaller  size  are  placed  in  region  one  near  the 
traction  free  edge  where  the  stress  distributions  are  of  the  greatest  interest. 


(a)  FOUR  LAYER  LAMINATE;  (0/90) g.  (90/0)3 


(b)  FOUR  LAYER  LAMINATE:  (45/-45)s 


(c)  EIGHT  LAYER  LAMINATE:  (90/0/-45/45)s,  (45/-45/0/90)s 


Fig.  4  Test  Cases 
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In  the  4-layer  cases  considered,  region  one,  where  =  .25  b,  is  subdivided 
into  25  equal  elements,  and  region  two,  where  J?.,  »  .75  b,  is  subdivided  into  5  equal 
elements;  This  mesh  was  chosen  on  the  basis  of  convergence  studies  which  illus¬ 
trated  that  the  use  of  additional  elements  showed  no  appreciable  change  in  the 
predicted  stresses  C93.  In  analyzing  the  8- layer  cases,  the  use  of  30  elements 
is  computationally  prohibative.  In  order  to  retain  accurate  stress  predictions 
near  the  free  edge,  the  size  of  region  one  is  reduced  such  that  2^^  »  .lb  and  » 
.9b.  Region  one  is  then  subdivided  into  5  equal  elements  and  region  two  is  sub¬ 
divided  into  6  equal  elements.  Additional  convergence  results,  based  on  a  4-layer 
laminate,  show  that  only  slight  changes  occur  in  the  stress  predictions  when  a 
more  refined  mesh  is  utilized.  No  subdivisions  are  needed  in  the  z  direction,  be¬ 
cause  the  element  is  developed  to  be  multilayered. 

The  total  width  of  the  ply  is  2b  while  the  height  of  each  layer  is  h.  The  ratio 
of  width  to  height  for  both  the  4 -layer  and  8 -layer  laminates  considered  is  4. 

The  boimdary  conditions  are  also  shown  in  figure  3b.  Symmetry  conditions  along 
y=b  (all  u  and  v*0)  and  along  z=0  (all  waQ)  are  imposed.  The  elastic  constants 
with  respect  to  the  principal  material  axes  of  each  ply  for  all  cases  considered 


are: 


=  20.0  X  10^  psi 
E22  =  2.1  X  10^  psi 
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“  '"si  *  ''^23 

G  - 

m  G- 

3  G 
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^31 
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0.21 

0.85  X  10^  psi 


The  first  case  considered  is  that  of  the  cross-ply  laminates,  CSO/OI^  and  CO/ 
903g' For  the  ctoes-ply  case,  the  u  displacement  is  a  function  of  x  only.  Conse¬ 
quently,  referring  to  the  expressions  for  the  v  and  w  displacement  given  in  (9), 
the  inplane  stresses  vanish  (i.e.  a  *  a  *  0).  Distributions  of  all  stresses 

A  Z 

(i.e.  (90°  layerh  (0°  layer),  cr^  (90°  layer),  (0°  layer),  c,  and  c^^) 

versus  y  along  the  0/90  interface  are  shown  in  figures  5  and  6  for  NO  SPECIAL  ele- 


Fig.  5  Stress  Results  for  Co/90]  laminate 


Fig.  S  (Continued) 
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Fig.  5  (Concluded] 
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Fig.  6  (Continued) 
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ment,  the  TFE  element  and  the  SC  element.  Also  shown  as  a  basis  of  comparison,  is 
a  solution  by  Wang  based  on  a  singularity  eigenfunction  expansion  with  botindary 
point  collocation  Cll]« 

For  both  the  CSO/O]^  and  [0/90]^  cases  all  three  approaches  agree  with  lamina¬ 
tion  theory  away  from  the  free  edge  which' is  expected.  The  differences  between  the 
solutions  become  apparent  as  the  free  edge  is  approached.  For  the  C0/903^  case 
(Figures  5a  through  5c) ,  the  solutions  for  in  the  90  and  0  layers  agree  for 
the  3  analyses.  The  solution  for  essentially  agrees,  although  for  the  TFE 
element  deviates  slightly  from  the  other  two  approaches  after  approximately  y/b  = 

.9.  Solutions  for  a  and  a  are  nearly  identical  for  the  NSE  and  the  SC  approaches, 
y  yz 

In  the  TFE  analysis,  these  stresses  are  forced  to  zero.  In  the  other  two  cases, 

both  a  in  the  90°  layer  and  a  tend  toward  zero  despite  the  fact  that  this 
y  y* 

condition  is  not  exactly  enforced.  The  value  of  in  the  0  layer,  however,  stays 
near  its  constant  value  near  the  free  edge  for  the  NSE  and  the  SC  approaches  and 
does  not  tend  toward  zero. 

Similar  results  are  observed  for  the  [90/0]^  cases  shown  in  figures  6a  through 

6c.  Again  the  results  for  are  similar  for  all  three  approaches,  and  the  results 

for  are  similar  although  slightly  different  after  approximately  y/b  =  .9.  Also, 

a  in  the  0°  layer  and  a  for  the  NSE  and  SC  element  again  tend  toward  zero,  al- 
y  y2 

though  they  are  not  forced  to  do  so  as  in  the  TFE  element  case,  whereas  in  the 
0°  layer  does  not. 

Notable  differences  between  the  [90/0]^  and  lO/90]^  cases  occur  in  the  results 
for  and  cr^^.  For  the  stacking  sequence  [0/90]^,  is  predominantly  negative 
and  rises  to  around  the  zero  point  near  the  free  edge,  while  for  the  C90/0]^  case, 
the  opposite  occurs;  however  the  two  curves  are  not  mirror  images,  for  the  Co/90jg 

case  dips  negative  between  about  y/b  «  .6  and  y/b  *  .9  and  then  rises  up  from  y/b  = 

.9  to  the  free  edge.  For  the  stacking  sequence  [90/01^,  J,  is  positive  up  to  about 
y/b  ■  .8  and  then  dips  negative,  rising  again  to  about  zero  at  the  free  edge. 
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The  important  differences  between  the  two  different  stacking  sequences,  and  the 
three  different  approaches  is  illustrated  by  looking  at  through  thickness  stress 
distributions,  a^,  and  are  plotted  through  the  thickness  in  Figures  7  and 
8.  Delamination  of  these  multilayer  structures  is  believed  to  be  caused  by  a  combi¬ 
nation  of  stresses  near  the  free  edge.  Consequently  these  distributions  are  of 
interest. 

For  the  C0/90j  laminate,  the  distribution  of  a  shows  about  the  same  shape  for 
s  z 

all  three  cases.  The  maximum  positive  value  is  larger  for  the  NSE  and  SC  analyses 

than  for  the  TFE  case,  where  as  the  maximum  negative  value  is  larger  for  the  SC 

and  TFE  approaches.  The  NSE  analysis  remains  negative  within  the  second  layer,  while 

the  other  two  cases  become  positive  and  then  become  zero  at  the  top  of  the  laminate 

along  with  the  NSE  approach.  The  stresses  a  and  a  shown  in  Figures  7b  and  7c 

y  yz 

respectively  are  forced  to  zero  in  the  TFE  element  case,  agrees  very  closely 
for  the  other  two  approaches,  tending  toward  a  large  negative  value  in  the  90°  lay¬ 
er  at  the  interface,  and  a  large  positive  value  in  the  0°  layer.  The  distribution 
of  oscillates  about  the  zero  point.  The  agreement  between  the  two  approaches, 

NSE  and  SC,  appears  poor,  but  note  that  is  an  order  of  magnitude  less  than  the 
other  stresses  in  question. 

The  C90/03  laminate  shows  a  much  different  distribution  of  a  .  Again,  however, 

S  2 

the  three  approaches  show  the  same  general  distribution  with  some  disagreement  in 
the  first  layer.  The  distribution,  however,  is  predominantly  compressive,  whereas 
the  distribution  of  for  the  C0/90jg  laminate  is  predominantly  tensile.  The  shape 
of  by  for  the  CSO/O]^  laminate  is  the  reverse  of  that  for  the  L0/90j^  laminate. 

The  two  appraches,  NSE  and  SC,  agree  well.  Recall  that  and  b^^  for  the  TFE  case 
are  set  equal  to  zero  along  the  free  edge.  The  distribution  of  b  again  oscillates 
about  zero,  and  is  an  order  of  magnitude  smaller  than  the  other  stresses. 

Recall  that  high  order  through  thickness  stress  distributions  have  been  assumed 

3  5  4 

within  each  layer;  b  is  of  order  z  ,  b  is  of  order  z  and  b  is  of  order  z  .  If 

y  •  z  yz 


Oz  VERSUS  z  AT  y»b  FOR  (0/90)s 

Through  thickness  stress  results  at  the  traction  free  edge  for  C0/90]  laminate 
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Fig.  7  (Continued) 
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(c)0  VERSUS  2  AT  y«=b  FOR  (0/90) 


(a)  O  2  VERSUS  z  AT  y»b  FOR  (QO/O)^ 

Fig.  8  ITjrough  thickness  stress  results  at  the  traction  free  edge  for  [90/0j  laminate 
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(b)a  VERSUS  z  AT  y-b  FOR  (90/0) 
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VERSUS  z  AT  y»b  FOR  (90/O) 
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lower  order  assumed  stress  distributions  could  accurately  predict  the  actual  stress 
distributions,  computation  time  could  be  reduced  without  effecting  the  results  of 
the  analysis.  If  it  is  possible  to  reduce  the  order  of  the  stress  assumptions,  from 
equilibrium  considerations,  each  stress  component  must  be  reduced  by  the  same  order. 
It  is  evident  from  Figures  7b  and  8b  that  the  order  of  can  not  be  reduced  with¬ 
out  loss  of  accuracy.  Therefore,  the  high  order  stress  distibutions  which  have 
been  used  are  in  fact  necessary. 

The  high  order  displacement  interpolations  might  also  be  lowered  to  save  on  com¬ 
putation  time.  Through  thickness  distributions  of  v  and  w  displacements  along  the 
free  edge  and  away  from  the  free  edge  for  the  SSE  and  SC  approaches  are  shown  in  Fig¬ 
ures  9a  through  9c.  In  this  analysis  the  v  displacement  is  of  order  and  the 

2 

w  displacement  is  of  order  z  .  The  v  displacement  is  chosen  one  order  higher  than 
the  w  displacement  so  that  their  re>ative  through  thickness  contributions  to 
will  be  of  the  same  order.  Away  from  the  traction  free  edge  (Figure  9c)  the  w  dis¬ 
placement  is  only  of  order  z.  To  be  consistent,  therefore,  the  v  displacement  would 
be  of  order  z^  although  it  appears  almost  constant  in  Figure  9c.  At  the  traction 
free  edge,  however,  it  appears  that  both  the  v  and  w  displacement  could  be  reduced 
by  one  order  of  z  without  appreciably  effecting  the  accuracy  of  the  displacement 
and  stress  distributions. 

The  next  case  considered  is  the  angle  ply  L45/-45^g  case.  In  this  case,  results 
by  Wang  and  Crossman  L4]  are  presented  as  a  basis  of  comparison  where  available.  In 
that  analysis,  constant  strain  triangles  were  used  to  model  the  problem.  Later  re¬ 
sults  by  Herakovich,  Nagarker,  and  O'Brien  C7]  ,  also  using  constant  strain  tri¬ 
angles,  show  results  which  differ  somewhat  from  Wang  and  Crossman  when  a  more 
refined  mesh  is  used.  The  results  for  this  case  are  presented  in  Figures  10a  through 

lOf.  Included  are  distributions  of  a,,  a  o  and  a  along  the  45/ -45  interface 

2  yz  xz  xy 

and  through  thickness  distributions  of  and  along  the  free  edge. 

Figure  10a  shows  the  distribution  of  a  along  the  45/-45  interface.  Results 
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(NO  SPECIAL  ELEMENT;  STRAW  CONTINUITY) 


Fig.  10  (Continued) 


(c)  ^VERSUS  y  ALONG  45/-45  INTERFACE  FOR  (45/-45) 
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(6)0^^  VERSUS  z  AT  y»b  FOR  (45/-46) 


ONO  SPECIAL  ELEMENT 

□tfe  element 


(0  O  2  VERSUS  z  AT  y-b  FOR  (45/-45)s 
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for  the  three  analyses  are  very  similar.  At  the  traction  free  edge  all  the  solutions 
reach  a  large  negative  value.  The  TFE  element  reaches  the  greatest  value  while  the 
NSE  and  SC  analyses  give  similar  results  at  the  traction  free  edge.  It  is  of  in¬ 
terest  to  note  that  the  Wang  and  Crossman  results  tend  to  a  large  positive  value 
while  all  of  the  results  presented  here  show  a  large  negative  value  at  the  traction 
free  edge. 


The  three  approaches  also  show  very  similar  results  for  the  distributions  of 

along  the  45/-45  interface.  Small  discrepancies  appear  near  the 

traction  free  edge  for  the  distribution  of  a  .  The  TFE  element  is  forced  to  have 

yz 

a  value  of  zero  at  this  point,  while  the  MSE  analysis  shows  a  small  positive  value 
and  the  SC  analysis  shows  a  small  negative  value.  The  results  for  show  that 
the  TFE  approach  again  gives  the  largest  value  at  the  free  edge,  while  the  results 
for  the  other  two  analyses  are  smaller  and  in  better  agreement.  Note  also  that  the 
results  for  the  SC  analysis  are  slightly  lower  than  the  other  approaches  from 
about  y/b  =  .6  to  the  traction  free  edge.  The  differences  observed  in  the  distribu¬ 
tion  of  0^^  between  the  three  approaches  is  also  near  the  traction  free  edge. 

Here  the  TFE  analysis  is  again  forced  to  be  zero.  It  also  shows  lower  values  from 
about  y/b  =  .35  to  the  traction  free  edge  when  compared  to  the  other  two  analysis. 
The  NSE  and  SC  analyses  give  essentially  the  same  results  but  do  not  reach  zero 


at  the  traction  free  edge. 


Th.e  three  approaches  again  show  similar  results  for  the  through  thickness  dis¬ 
tributions  of  a  and  u  at  the  free  edge.  The  values  of  a  through  the  thickness 
(Figure  10c)  are  all  positive.  The  TFE  element  reaches  a  greater  value  at  the  45/ 
-45  interface  than  either  the  NSE  or  SC  analyses,  but  the  relative  shapes  of  the 
distribution  for  the  three  appraches  is  the  same.  Again,  some  dissimilarities  are 
observed  between  the  three  analysis  for  the  distribution  of  o.  (Figure  lOf) .  They 
do  agree  well  at  the  midsurface,  where  the  ma.timum  positive  value  of  occurs.  The 
maximum  negative  value  occurs  at  the  45/ -45  interface,  where  the  TFE  element  again 


reaches  the  greatest  value  while  the  NSE  and  SC  approaches  tend  to  have  smaller 
negative  values. 


The  final  cases  analyzed  are  the  l45/-45/0/90j^  and  C90/0/-45/45j^  laminates. 
Again,  Wang  and  Crossman  results  are  included  on  the  graphs  where  available.  Fig¬ 
ures  11a  through  Ilf  show  the  results  for  the  r4S/-45/0/903  case.  Th.e  distribu- 

s 

tion  of  along  the  0/90  interface  and  along  the  midsurface  are  plotted  in  Figures 
11a  and  11b  respectively.  The  three  approaches  give  essentially  the  same  results 
up  to  the  traction  free  edge  where  the  NSE  analysis  drops  off  slightly.  The  results 
for  the  distributions  of  along  the  0/90  interface  (Figure  11c)  and  the  distri¬ 
bution  of  along  the  4S/-43  interface  (Figure  lid),  again  show  that  the  three 

analysis  agree  well  up  to  the  traction  free  edge.  The  distribution  of  shows 
the  NSE  approach  dropping  off  again  at  the  free  edge.  The  TFE  element  drops  down, 
but  then  is  forced  back  up  to  zero  at  the  traction  free  edge,  whereas  the  SC  approach 
rises  smoothly  to  zero.  All  the  results  for  rise  smoothly  at  the  free  edge.  The 
TFE  element  rises  to  the  highest  value  while  the  NSE  and  SC  approaches  rise  to  smal¬ 
ler  vail  e'  which  are  again  in  better  agreement. 

Through  thickness  distributions  of  cr  and  o  for  the  C45/ -45/0/90 j  laminate 
are  shown  in  Figures  lie  and  Ilf.  The  TFE  and  SC  approaches  show  tensile  stress  for 
7,  in  the  90°  and  0°  layers,  and  primarily  compressive  values  of  7^  in  the  -45°  and 
43°  layers.  The  NSE  analysis  gives  smaller  tensile  values  in  the  90°  layer,  but  re¬ 
mains  tensile  up  to  the  top  of  the  45°  layer.  Agreement  between  the  three  approaches 
is  slightly  better  for  a  .  Here  only  the  value  at  the  0/-45  interface,  where  the 
TFE  analysis  shows  a  larger  negative  value,  and  the  value  at  the  45/ -45  interface, 
where  the  TFE  analysis  shows  a  much  higher  positive  value,  differ  significantly  be¬ 


tween  the  three  analysis. 


Results  for  7^  along  the  0/90  interface,  and  7^  along  the  midsurface  are  shown  in 
Figures  12a  and  12b  for  the  C90/0/-4S/45j^  laminate.  Again  the  three  analysis  agree 
well  up  the  traction  free  edge.  The  NSE  approach  drops  off  to  give  a  negative  value 
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(a)  Oj  VERSUS  y  ALONG  0/90  INTERFACE  FOR  (45/-45/0/90) 
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Fig.  11  (Continued) 
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Fig.  11  (Continued) 


(e)  Oz  VERSUS  z  AT  y=b  FOR  ( 45/-45/0/90) 
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Stress  results  for  C90/0/-45/45n  laminate 


Fig.  12  (Continued) 
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for  along  the  0/90  interface  while  the  TFE  and  SC  analysis  rise  to  positive  val¬ 
ues.  Along  the  midsurface,  the  TFE  and  SC  approaches  reach  larger  negative  values 
than  the  NSE  anlaysis  at  the  free  edge. 

Distributions  of  a ^  and  through  the  thickness  at  the  traction  free  edge  for 
the  L90/0/-4S/4Slg  case  are  shown  in  Figures  12c  and  12d.  The  NSE  analysis  again 
shows  smaller  values  throughout  the  laminate  for  o  .  Just  as  in  the  l4S/-45/0/90j 

M  S 

case  it  does  not  show  the  more  definite  peaks  which  are  found  in  the  distributions 
for  the  TFE  and  SC  approaches.  The  results  for  the  TFE  and  SC  analysis  show  that 
is  compressive  throughout  the  first  three  layers  and  becomes  tensile  in  the  90° 
layer.  The  NSE  analysis  remains  compressive  up  to  the  top  of  the  90°  layer,  but  its 
value  remains  small.  The  distribution  for  a  demonstrates  the  same  behavior  as  that 

X2 

of  the  C4S/-4S/0/90]^  case  between  the  three  approaches.  Large  negative  values,  a- 

gain  for  which  the  TFE  analysis  shows  the  largest  occur  at  the  4S/-4S  interface. 

The  TFE  approach  also  shows  the  largest  positive  value  at  the  0/-43  interface  where 

all  three  approaches  reach  the  maximum  positive  value  through  the  thickness. 

The  differences  in  results  observed  between  the  L45/ -43/0/903  and  the  [90/0/ 

s 

-45/45 laminates  are  best  illustrated  by  considering  the  through  thickness  dis¬ 
tributions  of  o  _  and  a  .  The  laminate  stacking  sequence  where  the  90°  layer  is 
on  the  outside  shows-  predominantly  compressive  behavior  and  the  most  significant 
peaks  in  stress  are  compressive.  The  stacking  sequence  where  the  90°  layer  is  on 
the  inside,  however,  shows  predominently  tensile  stresses  and  all  significant  peaks 
are  tensile.  These  are  very  important  observations  considering  that  the  interlami¬ 
nar  stresses  near  the  free  edge  are  believed  to  be  the  cause  of  delaminatidn  in 
these  types  of  laminates. 


Ojxio'®ex 

(c)  O  2  VERSUS  2  AT  y»b  FOR  (90/0/-45/45) 
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(d)0^_  VERSUS  z  AT  y=b  FOR  (90/0/-45/45) 
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4.  Summary  and  Concluding  Remarks 

An  assumed  stress  hybrid  formulation  has  been  presented  for  the  problem  of  com¬ 
posite  laminates  under  \miform  inplane  strain.  A  special  purpose  multilayer  ele¬ 
ment  has  been  developed  which  satisfies  various  stress  and  strain  conditions  ex¬ 
actly.  The  multilayer  element  has  been  used  in  conjunction  with  three  basic  ap¬ 
proaches  to  the  problem;  Ij  NSE:  a  mesh  of  so  called  standard  elements  which  satis¬ 
fy  continuity  of  interlaminar  stresses  across  interlayer  boundaries  and  traction 
free  edge  conditions  along  the  top  and  bottom  of  the  laminate  are  used  thoughout 
the  mesh,  2)  7FE:  the  first  element  in  the  mesh  is  modified  to  satisfy  the  trac¬ 
tion  free  edge  condition  while  the  remainder  of  the  mesh  consists  of  the  standard 
element,  and  3)  SC:  the  standard  element  is  modified  to  satisfy  continuity  of  in¬ 
plane  strain  along  interlayer  boundaries  and  the  entire  mesh  is  comprised  of  these 
elements. 

Stress  results  for  fi\e  laminate  test  cases,  C90/o3  ,  Co/90j  ,  C45/-4S]  , 

S  5  S 

[l90/0/-45/45l]g  and  C4S/-4S/0/903g,  have  been  presented  and  discussed.  Basically 

the  three  approaches  show  consistent  results.  The  stress  contributions  which  are 

forced  to  tero  in  the  TFE  analysis  in  most  cases  also  tend  toward  zero  in  the  NSE 

and  SC  analysis  even  though  they  have  not  been  forced  to  zero.  Exceptions  to  this 

are  a  (0°  layer)  for  the  LO/90].,  C90/0J  laminates  and  a  for  the  1I45/-4S3 
/  ^  s  ^ 

laminate.  In  general,  the  TFE  analysis  displays  the  most  severe  distributions, 
exhibiting  the  highest  peaks  in  both  the  negative  and  positive  directions.  The  NSE 
approach  behaves  well  except  in  the  analysis  of  the  [90/0/45/ -45]^  and  [45/ -45/0/ 

90 jg  laminates  where  the  results  for  the  NSE  approach  drop  off  at  the  traction  free 
edge  and  therefore  do  not  agree  with  the  TFE  and  SC  approaches.  Looking  at  the 
through  thickness  plots  of  a,  and  a  along  the  traction  free  edge,  also  for  the 
[90/0/4S/-4Sng  and  [45/-4S/0/903g  laminates,  it  is  apparent  that  the  peaks  in  stress 
exhibited  by  the  TFE  and  SC  approaches  are  smoothed  out  by  the  NSE  anlaysis.  This 
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suggests  that  in  the  analysis  of  more  complicated  laminate  stacking  sequences, 
some  special  conditions  must  be  exactly  satisfied  in  order  to  obtain  con¬ 
clusive  stress  distributions. 

Some  general  observations  can  also  be  made  concerning  the  differences 
between  laminate  stacking  sequences.  It  is  observed  when  the  90°  layer  is 
on  the  outside  (i.e.,  both  for  [90/0]^  and  [90/0/45/-45] ^  laminates)  the 
through  thickness  distributions  for  and  (for  the  angle-ply  only) 
show  predominantly  compressive  stress  values.  Conversely,  when  the  90° 
layer  is  on  the  inside  the  distributions  for  these  stresses  are  predominant¬ 
ly  tensile.  This  suggests  that  a  laminate  stacking  sequence  with  a  90° 
layer  on  the  outside  will  be  less  likely  to  delaminate  under  an  inplane 
tensile  load  than  one  which  has  a  90°  layer  on  the  inside. 

Three  approaches  have  been  used  to  solve  the  problem  of  composite 
laminates  under  uniform  inplane  strain.  All  three  analysis  show  basically 
the  same  results.  Observed  differences  do  occur  in  the  vicinity  of  the 
free  edge.  It  is  not  possible  at  this  time  to  claim  that  one  or  the  other 
of  the  analysis  provides  the  correct  detailed  stress  distributions  for  the 
problem  in  question.  Conclusions  of  this  nature  must  await  an  independent 
analytical  solution  to  the  problem  which,  thus  far,  does  not  exist. 
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CHAPTER  3 

A  STUDY  OF  8 -NODE  SINGLE  LAYER  PLATE 
ELE>ENTS  WITH  A  STRAIGHT  TRACTION-FREE  EDGE 

ABSTRACT 

The  elements  developed  and  tested  in  the  last  chapter  were  based  on  a 
plane  2-D  theory.  For  general  multilayer  plate  problems  involving  free  edges, 
a  multilayer  plate  element  is  required  which  satisfies  the  traction-free 
conditions  along  one  of  its  edges.  The  present  chapter  describes  a  study 
which  is  needed  to  establish  the  basis  for  development  of  such  a  multilayer 
element.  Here  eight-node  single  layer  pure  bending  plate  elements  are 
developed  for  which  the  traction- free  conditions  are  exactly  satisfied 
along  one  straight  edge.  Transverse  shear  deformation  and  transverse  shear 
and  normal  stresses  are  included  so  that  the  elements  are  applicable  to  both 
thin  and  moderately-thick  single  layer  plates.  Various  plausible  stress 
fields  are  defined,  and  the  best  stress  field  (element)  is  identified  by 
comparison  of  results  for  selected  example  problems.  The  results  obtained 
in  the  computationally  efficient  pure  bending  study  can  then  guide  the  de¬ 
velopment  of  a  multilayer  element,  where  stresses  and  displacements  are 
assiimed  independent  within  each  layer.  In  such  an  element,  the  present 
stress/displacement  fields  must  be  extended  to  include  stretching  contri¬ 


butions. 


1.  Introduction 


In  the  analysis  of  multi-layer,  laminated,  composite  plate  structures,  it  has 
been  observed  [1,2]  that  near  the  traction- free  edges  of  such  structures  severe 
gradients  in  the  interlaminar  stresses  may  exist.  These  gradients  can  lead  to  de¬ 
lamination  as  well  as  other  forms  of  laminate  failure.  Therefore,  it  is  necessary 
that  accurate  representation  of  the  stress  fields  near  traction-free  edges  be  ob¬ 
tained.  In  employing  the  finite  element  method  to  solve  this  class  of  problems,  ac¬ 
curate  analysis  in  the  vicinity  of  traction-free  edges  appears  to  require  the  de- 
velopement  of  a  special  purpose,  multi-layer,  plate  element  which  exactly  satisfies 
the  traction-free  conditions  along  at  least  one  edge.  (i.e.  normal  and  shear  stresses 
are  sero  along  one  edge  of  the  element.) 

Historically,  plate  elements  have  been  based  on  the  assumed-displacement  formula¬ 
tion.  In  this  approach,  displacement  boundary  conditions  are  exactly  satisfied  while 
stress  conditions  are  satisfied  only  in  an  approximate  (weighted,  integral)  sense. 
Alternatively,  in  the  assumed-stress  hybrid  formulation  both  stress  and  displace¬ 
ment  boundary  conditions  can  be  satisfied  exactly.  The  hybrid-stress  model  is  a  two- 
field  principle  in  which  equilibriating  intraelement  stress  fields  and  compatible 
displacement  fields  are  assumed  independently.  The  stress  parameters  are  eliminated 
on  the  element  level  and  a  conventional  stiffness  matrix  results  [3].  In  view  of 
the  independent  assumption  of  stresses  within  each  element,  it  is  possible  to  ex¬ 
actly  satisfy  traction-free  conditions  by  appropriate  choice  of  the  stresse  fields. 
This  feature  establishes  the  assumed-stress  hybrid  formulation  as  a  viable  approach 
for  developing  such  special  purpose  elements.  In  many  cases  [4,S],  ,hybrid-stress 
elements  have  been  found  to  yield  improved  convergence  and  intraelement  stress  pre¬ 
dictions  in  comparison  with  analagous  assumed-displacment  elements. 

The  hybrid-stress  model  is  also  well  suited  to  the  development  of  multi-layer, 
laminated,  composite,  plate  elements.  Stresses  and  displacements  may  be  assumed  in- 
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dependently  within  each  layer  and  appropriate  interlayer  displacment  and  stress 
continuity  conditions  satisfied  exactly  C6,7l.  Such  elements  have  been  shown  to  be 
applicable  for  up  to  moderately  thick  laminates  Transverse  shear  effects  can 
also  be  included  in  a  less  general  manner  (i.e.  on  the  laminate  level  CqD)*  These 
hybrid-stress  elements  have  been  found  to  be  more  accurate  than  comparable  assumed- 
displacment  elements  Cl0,ll3. 

Based  on  these  observations,  the  hybrid-stress  model  appears  to  be  the  ideal  choice 
for  the  development  of  a  special  purpose,  multi-layer,  traction-free  edge,  plate 
element.  However,  before  this  element  can  be  developed,  it  is  necessary  to  examine 
the  possible  stress  fields  assumed  in  the  element  interior  which  satisfy  the  trac¬ 
tion-free  conditions  along  an  edge.  This  is  best  accomplished  by  first  considering 
single-layer,  isotropic,  plate  elements.  Once  the  best  stress  fields  are  identified 
for  single-layer  plates;  the  results  can  then  be  extended  to  multi-layer  plates. 

In  general,  the  elements  cited  earlier  in  this  discussion  have  been  based  on  4- 
node  bilinear  displacement  fields.  But,  in  a  recent  series  of  articles  Cl2,13,14]] 
a  family  of  single-layer,  isotropic,  plate  elements  based  on  the  hybrid-stress  mod¬ 
el  have  been  developed  and  tested.  These  elements  use  displacement  distributions 
based  on  Mindlin  plate  theory  ClSH  and  include  all  components  of  stress.  Their  ad¬ 
vantage  over  analagous  assumed-displacement  elements  is  that  the  element  stiffness 
matrix  exhibits  correct  rank  and  accurate  solutions  can  be  obtained  for  arbitrarily 
thin  to  moderately  thick  plates.  That  is,  difficulties  regarding  -excessive  stiffen¬ 
ing  (i.e.  locking)  observed  in  Mindlin-type  assumed-displacment  elements  are  avoided 
in  the  hybrid-stress  elements.  Similar  advantages  should  be  expected  if  multi-layer 
versions  of  these  elements  are  developed  for  the  analysis  of  laminated,  composite 
plates. 

In  the  present  study  a  single-layer,  isotropic  plate  element  is  developed  for  which 
traction-free  conditions  are  exactly  satisfied  along  one  edge.  After  choosing  an  ap¬ 
propriate  displacement  field,  various  plausible  intraelement  stress  fields  are  de- 
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fined.  The  resulting  special  purpose  elements  are  coupled  with  the  compatible  stan¬ 
dard  element  to  perform  the  analysis  of  several  plate  problems  which  include  a  trac¬ 
tion-free  edge.  The  performance  of  the  special  purpose  elements  is  evaluated  and 
their  potential  value  is  assessed  for  the  development  of  a  special  purpose,  multi¬ 
layer,  plate  element. 


Stress  Interpolations 


A.  Basic  Equations 

To  establish  a  basis  for  the  present  work,  the  assumed-stress  hybrid  formulation 
for  plate  C12D  will  be  summarized  here.  The  hybrid-stress  functional  can  be  writ¬ 


ten  as: 


w  I  j  %  I  -  I  o^'^dV  ♦  I  u'^T  ds 

n  *n 


where: 


stress  vector 

strain  vector  calculated  from  displacements,  u. 

displacement  vector 

prescribed  tractions 

material  property  matrix 

volume  of  the  nth  element 

boundary  of  the  nth  element  over  which  tractions 
prescribed 


Based  on  Mindlin  plate  theory  Cl53,  the  through-thickness  displacement  distri¬ 
butions  are  assumed  in  the  form  (pure  bending  only) : 


u(x,y,z)  •  rey(x,y)  ^2) 

v(x,y,z)  -  -rejj(x,y) 

w(x,y,z)  •  w(x,y) 

-where  positive  sign  conventions  for  displacements  and  rotations  are  shown  in  Figure 

1. 

From  (2),  the  generalized  displacements  0  (x,y),  0  (x,y),  and  w(x,y)  can  be  ex- 
pressed  in  terms  of  a  set  of  nodal  degrees  of  freedom  0^^^,  0^^,  and  w^  by  construc¬ 
ting  a  set  of  continuous  shape  functions  to  use  as  displacement  interpolations. 
Applying  the  linear  strain-displacement  relations  yields: 
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where  q  is  the  vector  of  nodal  degrees  of  freedom  for  an  element,  and  the  sub¬ 
scripts  "f"  and  "s"  refer,  respectively,  to  flexural  and  transverse  shear  con¬ 
tributions. 

The  stresses  are  expressed  polynomial  form  in  terms  of  a  set  of  stress  para¬ 
meters,  6.  The  stress  assiauption  is  required  to  satisfy  the  3-D  homogenious 
equilibrium  equations.  Moreover,  for  a  plate  loaded  transversely  at  :=h  (see 
Figure  13.the  free-surface  conditions  can  be  expressed  ass 

-  0 

Oyj(x.y>h)  •  0  (4) 

o^Cx.y^-h)  ■  0 


Based  on  the  equilibrium  equations  and  the  free-surface  conditions  of  (i) ,  the 
through-thickness  distribution  of  stresses  is  assumed  in  the  form  (corresponding 
to  pure  bending  contributions  only) s 


Froa  (S),  it  Is  observed  thst  first  the  polynaalsls  for  9^,  <7^  ere  defined. 
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after  which  c,,  5"^,  a  are  calculated  froo  the  last  half  of  (S).  Note  that  this 
x^i  yx  z 

gtiarantees  satisfaction  of  the  3>D  homogeneous  equilibrium  equations.  The  stress 
assumption  can  then  be  cast  in  the  form: 


Substituting  (3)  and  (6)  into  (1)  and  manipulating  the  result  (described  in  de¬ 
tail  in  C12D  )  yields  the  stiffness  matrix  for  the  plate  element: 


where: 


i  -  C%-)  CGf  ♦  -  %-5.5‘^CGf  ♦ 

■[  !f  5f  "[  li  5 

S.  *  1  ?*  5> 


dA 


dA 

r  1 

-w 


S  A. 


If  -  1/2 


-V 


V  Zh‘ 

u  ? 


0 

0 


V  2h^  V  2h^ 


S2  h  0 

13? 


ZCl-^v)  J, 


(7) 


(8) 


C9D 


and  A.  is  the  area  of  the  midsurface  of  the  nth  element, 
n 

To  implement  this  nuterically,  it  is  common  practice  to  combine  the  flexural  and 
transverse  shear  parts  such  that:  G  ■  G£ 

H  .  ^Th£  +  —  H  ) 

-  3  S 


(10) 


where? 


G  -  J  B 

H.v  f*!  Jd 


n. 


Then,  the  stiffness  matrix  is  simply: 


.  k  -  gTh"^g 


Clla) 

(lib) 


(12) 


The  integrals  in  (11)  are  mapped  from  the  x-y  plane  into  the  C-n  plane  and  nimier- 
ically  integrated  by  Gauss  quadrature. 

With  this  information  as  a  basis,  we  will  next  look  into  the  formulation  of  the 
displacement  interpolation  (4  ■  Bq)  and  the  stress  assumption  (?  •  FS)  necessary 
to  yield  a  single-layw,  traction-free  edge,  plate  element. 

B.  Displacement'  Interpolation 

In  the  work  presented  in  references  Cl2,13,14!],  a  family  of  isoparametric  plate 
elements  based  on  the  assumed-stress  hybrid  formulation  and  using  the  Serendipity 
family  of  displacement  interpolations  was  developed  and  tested.  Comparison  of  the 
elements  in  the  family  shows  the  8-node  and  12-node  elements  to  be  more  accurate 
per  degree  of  freedom  than  the  4-node  element  Cl4n.  However,  12-node  elements  are 
generally  perceived  as  too  complex  for  practical  applications  and  in  applications 
such  as  nonlinear  analysis  where  computation  time  is  strongly  dependent  on  element- 
level  operations.  In  light  of  this,  the  8-node  element  (element  QHl  Cl3D  )  was  cho¬ 
sen  to  interface  with  the  special  purpose  element. 

In  order  to  insure  that  displacements  are  continuous  between  elements,  the  dis¬ 
placement  interpolation  of  the  special  purpose  element  must  be  identical  to  the 
displacement  interpolation  of  the  8-node  element;  that  is,  the  quadratic,  Serendip¬ 
ity  shape  functions.  In  the  C-h  plane,  they  are  in  the  form: 
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e_C5.n3  -  I  N. cc.n)e 

! - ^ - Ji- 

e^C5.n)  -  I  H,(5.n30 

y _ i«i  ^  _ _ 7^ 

where  are  the  shape  functions  and  3^,  6^^,  are  the  degrees  of  freedom  (gen¬ 
eralized  displacements]  At  node  i. 

For  isoparaaetric  elements,  the  coordinate  mapping  is  of  the  same  form  as  (13] 
(i.e.  quadratic];  but,  for  the  special  purpose  element,  only  bilinear  mapping  is 
allowed.  This  permits  it  to  take  on  a  general,  quadrilateral  shape,  but  insures  that 
its  sides  remain  straight.  (This  is  a  requirement  of  the  stress  assumption,  and  will 
be  explained  shortly].  So,  the  mapping  is  given  by  the  bilinear  shape  functions: 

4” 

I  Ni(5»n]x. 

'  i«l  ^  ^ 

■  (14] 

y  •  I  i^iCC.nJy, 

i-i  -  - 

where  (x^,  y^^]  are  the  coordinates  of  node  i. 

Recasting  (3]  into  the  plane  results  in: 


f^x  3  3x  3  \ 
SrT  ”  5^  S?"/ 


(iL  3 


in) 


•  3^  3v  3  ^ 

ySrT  5^  5^  5^ ! 


^3x  3  3x  3  \ 

Srf  Sn"  j 


^  Ini 


'  w(5.n] 
9^(5  .n] 

9y(5.n] 


3  3x  3  \ 

5n  *  5n  JT; 


36 


where  |-j|  is  the  Jacobian  of  the  coordinate  transformation. 

Substituting  (13)  into  (IS)  yields  the  B  matrix  (fr«  Bq)  which  is  used  in  the  in¬ 
tegral  of  (11a). 

C.  Stress  Assumption 

The  choice  of  a  stress  assumption  for  the  special  purpose  element  is  not  obvious. 
In  this  section,  several  plausible  stress  fields  are  developed.  Subsequently,  they 
will  be  subject  to  numerical  tests  in  order  to  identify  the  better  stress  fields. 

First,  consider  a  rectangular  elaaent  with  sides  parallel  to  the  x  and  y  axes  as 
shown  in  Figure  2a.  Letting  x>0  be  the  traction-free  edge  leads  to  the  requirement 
that: 

0^(0, 7.2)  -  0 

a^(0,y,z)  «  0  (16) 

Since  these  conditions  must  hold  for  all  values  of  s,  the  relations  in  (5)  can  be 
used  to  specifically  require  that: 

--  (I7a). 

°  Cl7b) 

axjCO.'y)  ^  C17c) 

Notice,  that  these  conditions  simply  require  that  a  ,  a  ,  and  o  be  zero  along 

A  'O' 

the  line  x>0  for  all  y.  For  all  other  x  and  y  the  stresses  can  vary  according  to 
the  polynomials  vdiich  represent  them  in  the  stress  assumption.  In  fact,  for  an  ele¬ 
ment  of  general,  quadrilateral  shape,  it  is  possible,  in  a  finite-element  analysis, 
to  define  a  local  coordinate  system  for  that  element  such  that  its  traction-free 
edge  is  on  the  line  x>0  (see  Figure  2b).  Then  the  conditions  of  (17)  can  be  applied 
(replacing  x  by  x  and  y  by  y)  in  the  local  system  and  k  can  be  formed  in  that  sys¬ 
tem.  Finally,  k  can  be  transformed  into  the  global  system  by  the  standard  transfor¬ 
mation  of  displacements  [163.  Furtheimore,  equation  (17)  requires  0^^,  and 
to  be  zero  on  the  line  x  >  0  (x-0)  which  is  a  straight  line.  Therefore,  the  sides  of 


(B) 


Figure  2.  Typical  Traction-Free  Edge  Plate  Elements 
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the  element  must  remain  straight  in  order  for  these  conditions  to  apply.  This  ex¬ 
plains  why  the  coordinate  moping  was  defined  to  be  bilinear:  to  allow  the  element 
to  take  on  a  general,  quadrilateral  shape,  but  insure  that  its  sides  remain  straight. 

Another  criterion  which  the  special  purpose  element  is  subject  to  is:  the  number 
of  stress  parameters,  n^,  must  be  greater  than  or  equal  to  the  number  of  degrees 
of  freedom,  n^,  minus  the  number  of  rigid  body  modes,  n^,  i.e. : 


"S  >  "q  ’  "r 


(18) 


This  relation  is  a  necessary  but  not  sufficient  codition  to  guarantee  an  element 
stiffness  matrix  of  correct  rank  [Il73. 

To  summarize,  the  stress  assumption  for  the  special  purpose  element  is  defined  by 

first  defining  the  interpolations  for  a  ,  a  ,  and  a  as  functions  of  in-plane  co- 

X  y  xy 

ordinates,  x  and  y,  such  that  the  conditions  specified  by  (17)  and  (18)  are  met.  The 

forms  of  the  remaining  stress  components,  o  ,  a  ,  and  o  ,  are  then  determined  frcm 

XZ  /  z  z 

the  last  three  equations  of  (5).  This  guarantees  that  the  equilibrium  equations  aiid 
the  free-surface  conditions  of  (4)  are  satisfied  exactly. 

As  a  starting  poirjt,  consider  the  following  stress  field: _  ^ 

“  h  *  hoy^  * 


"xy  “  ®16  *  ^2Z^^y  * 

827XV  +  *  330^^ 


(19) 


"  ®31  ^  *  ^ZZy  *  ^ZA^y  *  ^ZS^^  *  ^36^^  *  ^37^^^  ^  ^ZZ^^y 


Note  that  a  and  0  are  full  qnartlc  polynomials,  while  a  is  identical  to  that 
X  X/  y 

used  in  element  QHl  Cl^D  (a  reasonable  starting  point  since  is  not  constrained 
by  (17)). 

Applying  equations  (17a),  (17b),  and  (17c)  (after  calculation  of  a_.)  yields,  re- 

aZ 
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spectively: 


61-83-85-810-815  -  0  .• 

(20a) 

*l6'*l8'*2l'hs’^0  *  . 

(20b) 

,  ■  - -  — 

(20c) 

®2*®S"®9"®14  "  ° 

After  imposing  equations  (20)  and  renumbering  the  the  stress  field  becomes: 

34*^  ♦  65x^7. 

V"  ®8*y  "  Sgxy^  *  B^^xy^  .  B^jX^  *  *  B^jxV  .  Si4X^.  S^^y  > 

V  “  ®7  ■"  ®8y  "■  *  ®loy^  *  ^BjjX  ♦  28^2*/  ♦  28^3X7^  ♦  38^^x^  ♦ 

46i6X^  *  8i9  *  620*  *  2822y  i' 2823x7  *  S^x^"'" 


c,  «^2Sj  ♦  262/  •*•  2837^  +  66 4X  ♦  665X7  ♦  128 gX^  ♦  26g  *  48^7  *  SSj^y^  +  '^®12*  *' 
88^5X7  ♦  66j5X^  •*•  2622  ^®23* 


’y  ■  *17  *  *18*  *20*>'  *  *21*^  *  *22>'^*  *23*>'^*  *24*^>' 


*  282x7  ♦  283x7*^  ♦  SB^x'^  +  SBgX'^y  +  46gX^  +  8gX  ♦  28gxy  +  SB^^xy 

*12*^  •  **13**>'  •  *13*’ 


J.- 


-  The  3-q  relation  of  equation  (18)  for  the  8-node  plate  elements  is: 

"a  >  n-  -  «,  ■  24  -  3  ■  21  (22) 

p  •  r 

Therefore,  with  24  stress  parameters,  the  stress  assumption  given  by  (21)  satis¬ 
fies  the  criteria  established  in  equation  (17)  and  (18)  for  the  special  purpose  ele¬ 
ment.  Equation  (21)  can  be  cast  in  the  foxn  of  (6)  to  yield  the  7  matrix  which  is 
used  in  the  integrals  of  (11). 

The  single-layer,  traction-free  edge  plate  element  based  on  the  24-8  field  of  (21) 
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will  be  denoted  as  element  PL24. 

In  any  finite-element  formulation,  a  necessary  first  test  of  an  element  is  to 
check  its  rank  by  calculating  the  eigenvalues  of  the  element  stiffness  matrix,  k. 

The  number  of  zero  eigenvalues  should  equal  the  number  of  rigid  body  modes  in  the 
element.  Any  additional  zero  eigenvalues  correspond  to  spurious  zero  energy  modes 
(i.e.  additional  kinematic  modes,  AXM)  which  must  be  eliminated  or  constrained  be¬ 
fore  the  element  can  be  safely  used  in  a  general  finite-element  analysis  ClSD. 

Eigenanalysis  of  k  for  element  PL24  shows  S  zero  eigenvalues.  With  3  rigid  body 
nodes,  this  indicates  the  presence  of  2  AKM.  Inspection  of  k  reveals  that  the  ro¬ 
tational  degrees  of  freedom  at  the  center  node  along  the  traction-free  edge  (i.e. 
node  8  in  Figure  2)  receive  no  stiffness  contributions.  This  result  is  more  appar¬ 
ent  by  considering  the  more  conventional  expression  for  G  (equation  (11a))  calcu¬ 
lated  as  an  integral  of  tractions  (from  stresses)  times  displacements  along  the  ele¬ 
ment  surfaces.  All  tractions  (traction-free  edge  and  upper/lower  surfaces)  which 

multiply  9  and  9  (rotational  degrees  of  freedom  at  node  8)  are  zero.  Hence,  the 
y  ®  xs 

columns  in  G  corresponding  to  9  and  9  are  zero,  and  the  corresponding  rows/col- 
••  y  8  X0 

umns  in  k  will  be  zero.  Since  these  2  degrees  of  freedom  have  no  stiffness  associ¬ 
ated  with  them,  they  can  be  eliminated  from  k.  Numerically,  this  is  done  by  artifi¬ 
cially  constraining  these  degrees  of  freedom  to  be  zero.  Note  that  this  operation 
is  equivalent  to  redefining  the  interpolations  for  9  and  9  to  be  linear  along  the 

/  A 

traction-free  edge.  This  resulting  stiffness  matrix  essentially  has  two  fewer  de¬ 
grees  of  freedom  so  that,  in  this  case,  n^  »  24  -  2  ■  22.  This  changes  the  require¬ 
ment  on  number  of  stress  parameters  given  by  (22)  to: 


no  >  n^  -  n  ■  22  -  3  ■  19 

B  q  r 


(23) 


A  subsequent  eigenanalysis  of  k  for  element  PL24  with  9  and  9  constrained  yielded 

^  /  8  a8 

3  zero  eigenvalues  corresponding  to  the  3  rigid  body  modes.  With  the  2  AKM  elimi¬ 
nated  through  the  constraint  of  Sy^and  9^^,  element  PL24  is  a  viable  candidate  for 
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use  as  a  single-layer,  traction-free  edge,  plate  element. 

Numerical  experience  of  several  authors  Cl9,20D  suggests  that  the  "optimum"  num¬ 
ber  of  6.  is  near  the  minimum.  Since  n^  ^  19  and  n.  ^  24  in  element  PL24,  it  seems 

X  DP 

reasonable  to  explore  ways  to  reduce  the  number  of  8^^.  Eliminating  certain  8^^  can 
result  in  the  introduction  of  AXM,  however.  To  deteimine  which  8^^  can  be  safely  elim¬ 
inated,  one  must  investigate  the  equation: 

G  a  =  0  C24) 


where  G  is  given  in  (11a)  and  a  is  a  vector  of  generalized  displacement  parameters 
(coefficients  of  the  polyncxnial  interpolations  for  displacement  which  can  be  unique¬ 
ly  related  to  the  actual  degrees  of  freedom)  for  an  element.  The  solution  a  «  0 
corresponds  to  the  rigid  body  modes;  any  other  non-trivial  solution  corresponds  to 
AKM  [21]. 


Equation  (24)  was  evaluated  for  element  PL24  using  a  square  of  side  length  2  (this 
is  a  sufficiently  critical  geometry  for  investigating  AXM) .  Based  on  the  result, 
the  following  observations  were  made: 


1)  8g  and  8]^q  are_redundant ;  can  be  eliminated  since  it  represents  a  higher 
order  term  in  a  . 


ar- 
order 


2)  2  of  the  following  6.  can  be  eliminated:  8.,  8,,  8^,  8,,  6-, 
guments  of  completeness,  only_consider  eliminating^two  of  the  highest  ord 
terms:  8j,  8^  in  °x*  ®15  in  o^. 

3)  2  of  the  following  8.  can  be  eliminated:  8g,  8,,,  8,^,  8,,.  Again^^  only  consi¬ 
der  eliminating  two  of  the  highest  order  termsT'^Bj^j,  8j^^,  8j^g  in  o  . 


xy 

With  no  further  information,  the  task  of  eliminating  6^  from  (21)  and  identifying 
the  "best"  traction-free  eleement  is  one  of  trial  and  error.  Numerous  candidate 


stress  fields  can  be  defined  by  elmination  of  combinations  of  stress  parameters 
(while  preserving  correct  stiffness  rank) ;  however,  certain  of  these  fields  have 
been  eliminated  on  the  basis  of  preliminary  numerical  experimentation.  The  observa¬ 
tions  made  in  these  preliminary  tests  will  be  briefly  summarized. 

First,  from  1),  setting  ■  0  in  (21)  yields  a  23-8  element  which  behaves  al¬ 
most  identically  to  element  PL24.  This  is  eiipected  since  8^g  represents  a  redundant 
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equation  in  the  solution  of  equation  (24) .  (For  this  reason,  will  not  be  con¬ 
sidered  in  any  subsequent  elements.)  Second,  removing  all  the  quartic  terms  from 
(i.e.  ®io"®13*®15"®16  "  results  in  a  20-6  element.  This  element  has  correct 
rank  as  expected;  but  yielded  poor  moment  distributions  in  the  preliminary  tests. 

It  therefore  seems  advisable  to  have  some  quartic  teims  in  .  Based  on  these 

xy 

observations,  the  candidate  stress  asstanptions  are  reduced  to  the  following  few 
which  will  be  described  and  assessed  in  more  detail: 

4  — 

ELEMENT  PL21:  Remove  the  highest  order  terms  present  in  equation  (21)  (x  in 
and  a^y)  by  setting  *  0.  (Also,  6jq  *  0) 

ELEMENT  PL19:  Since  the  traction-free  edge  condition  of  (17)  removes  the  constant 
and  linear  terms  in  x  from  and  only  the  constant  terms  in  x  from 
a^y,  it  seems  reasonable  to  expect  to  be  one  order  higher  in  x  than 
a  .  Therefore,  by  removing  B,  from  a  to  make  it  cubic  in  x,  c 

o  X 

should  be  reduced  to  a  quadratic  form  in  x  by  removing  B^^,  B^^,  and 
Bjg.  So,  set  86-6i4*8i5=Bjg  »  0  and  Bjq 

ELEMENT  PL20:  The  traction  free-edge  condition  forces  —  as  well  as  to  be  zero 
along  the  free  edge.  With  these  strict  constraints,  may  have  to  be 
of  thehighest  order  possible  in  x  in  order  to  acciirately  predict  the 
stress  field.  So,  referring  to  equation  (21),  retain  all  terms  in 
and  remove  Bj^g,  and  B^g  from  a^y  as  in  elmnent  PL19;  i.e.:  6j^^= 

®15"®16"  °  ®10  “  °- 

Though  several  other  stress  assuonptions  were  considered,  the  elements  described 
above  are  the  best  candidates  to  use  as  single-layer,  traction-free  edge  plate  ele¬ 
ments. 

In  closing,  it  should  be  mentioned  that  attempts  were  made  to  develop  traction- 
free  elements  with  arbitrary  curved  traction-free  edges.  In  this  approach,  the  ele¬ 
ments  were  mapped  into  ^-ri  space  where  the  stress  assumptions  were  to  be  applied. 
However,  the  stress  assumptions  are  required  to  satisfy  the  3-D  equilibrium  equa- 
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tions;  vrtiich,  in  this  case,  had  to  be  written  in  a  non-orthogonal  system,  since  the 
elements  had  arbitrarily  curved  edges.  Using  a  tensor  analysis  approach,  it  was 
found  that  writing  a  stress  assumption  which  satisfied  the  equilibrium  equations  in 
a  non-orthogonal  space  was  not  only  an  extremely  long  and  difficult  task,  but,  could 
not  guarantee  greatly  improved  results  for  the  additional  effort.  In  fact,  as  will 
be  seen  in  the  next  chapter,  approximating  ctirved  edges  with  straight  segments  yields 
good  results  in  most  cases. 
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3.  Example  Problems  and  Numerical  Results 

In  the  last  section,  several  plausible  single-layer,  traction- free  edge,  plate 
elements  were  developed. (Summarized  in  Table  1).  All  the  elements  satisfy  the  trac¬ 
tion-free  edge  conditions  of  (17)  as  well  as  the  3-q  relation  of  (18).  Also,  the 
stiffness  matrix  of  each  element,  k,  is  of  correct  rank,  (i.e.,  three  zero  eigen¬ 
values  corresponding  to  rigid  body  motion  after  elimination  of  the  two  AKM  associated 

with  9  and  0  ) . 

ys  xa 

In  order  to  determine  the  better  elements,  it  is  necessary  to  assess  their  per¬ 
formance  in  the  numerical  solution  of  a  few,  selected,  example  problems.  This  will 
not  only  identify  the  better  special  purpose  elements,  but,  determine  whether  a  need 
exists  for  such  elements  in  the  analysis  of  single-layer,  isotropic  plates.  More¬ 
over,  this  may  provide  insight  as  to  the  potential  value  of  these  elements  in  subse¬ 
quent  development  of  multi-layer,  traction-free  edge,  plate  elements. 

In  this  section,  two  example  problems  are  considered.  The  performance  of  the 
special  purpose  elements  of  Table  1  is  compared  with  exact  solutions  as  well  as  with 
the  non-traction-free  edge  element,  QHl.  (Note,  that  element  PL24  is  not  considered 
since,  as  mentioned  earlier,  it  has  too  many  stress  parameters) . 

The  first  example  problem  is  shown  in  Figure  3.  It  is  a  square  plate  (length, 

L;  thickness,  2h)  with  two  opposite  edges  simply-supported  and  the  other  two  free, 
subject  to  a  transverse  uniform  load,  p[21].  By  symmetry,  it  is  possible  to  model  a 
quarter  of  the  plate  in  the  FEM  analysis  using  an  NxN  mesh  of  elements.  The  special 
purpose  elements  are  situated  along  line  AC.  Results  are  presented  for  the  case: 

L=10.0  in,  h=0.05  in,  p^=5.0  psi,  E»3.0xl0^  psi,  and  v*0.3. 

In  Table  2  deflections  at  the  center  and  left  edge  (free  edge)  of  the  plate  are 
compared.  Notice,  that  even  for  the  coarsest  mesh,  deflections  are  within  4%  of  the 
exact  solution,  and  as  the  mesh  is  refined  the  deflections  converge  to  less  than  1%  erroi 
Also,  in  the  coarser  meshes,  the  special  purpose  elements  predict  displacements  slight- 


ELEMENT 


ZEROED  P’S 


NUMBER  OF  p'S 


PL19  p6,p^0,p^4,p^5,p^e  19 

PL20  ^10,^14,pi5,i?16  20 

PL21  p6,P^0,P^e  21 

PL24  24 


NOTE:  P’S  BASED  ON  EQUATIONS  (21) 


Table  1.  Single-Layer,  Traction-Free  Edge,  Plate  Elements 


UNIFORM  LOAD.  P 


Figure  3.  Problem  1 .  Rectangular  Plate 
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1x1  MESH 


El£MENT 

Wq 

TERROR 

Wc 

EXACT 

EEM 

EXACT 

fEM 

-0  fennUn 

QH1 

0.23824 

0.2407 

•1.033 

0.27464 

0.2848 

-3.699 

PL19 

0.2403 

-0.865 

0.2826 

-2.898 

PL20 

0.2372 

0.437 

0.2753 

•0.240 

PL21 

1 

t 

0.2373 

0.395 

1 

0.2753 

-0.240 

2x2  MESH 


ELEMENT 

Wo 

TERROR 

Wc 

EXACT 

FBM 

EXACT 

REM 

%  ERROR 

QH1 

0.23824 

0.2386 

-0.151 

0.27464 

0.2776 

•1.078 

PL19 

0.2384 

-0.067 

0.2772 

•0.932 

PL20 

0.2383 

-0.025 

0.2750 

-0.131 

PL21 

1 

f 

0.2383 

-0.025 

f 

0.2750 

•0.131 

4x4  MESH 


ELEMENT 

Wq 

•b  ERROR 

Wc 

"oERROfl 

EXACT 

FEM 

EXACT 

FEM 

QH1 

0.23824 

0.2384 

-0.067 

0.27464 

0.2753 

-0.240 

PL19 

0.2384 

-0.067 

0.2753 

-0.240 

PL20 

0.2384 

-0.067 

0.2743 

0.124 

PL21 

1 

r 

r 

NOTE:  %ERR0Bs(1-FEM/EXACT)X  100% 


Table  2.  Problem  1.  Deflection  Comparisons 
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ly  better  than  element  QHl;  but,  in  general,  all  of  the  elements  converge  quite 
well  and  the  special  purpose  elements  show  no  particular  advantage  at  this  point. 

The  largest  moments  and  stresses  in  Problem  1  occur  along  line  CD.  (Note:  See 
Figure  1  for  sign  convention  for  moments  and  stresses).  To  further  assess  the  ele¬ 
ments,  Problem  1  was  run  with  a  2x2  mesh  and  the  moments,  M  and  M  ,  along  line  CD 

y  X 

were  plotted  (Graphs  1  and  2,  respectively).  along  CD  is  a  maximum  at  the  free 
edge  and  decreases  to  a  normalized  value  at  the  plate  center.  Notice,  that  all  the 
special  purpose  elements  show  some  deficiency:  either  they  over  estimate  the  moment 
at  the  free  edge  or  underestimate  the  moment  at  the  interface  of  the  special  pur¬ 
pose  elements  and  regular  elements.Of  the  special  purpose  elements,  PL19  provides 
the  best  solution;  however,  the  standard  element,  QHl,  also  yields  a  reasonable  pre¬ 
diction. 

is  calculated  from  and  so  is  zero  at  the  free  edge.  Along  line  CD,  it  goes 
from  zero  at  the  free  edge  to  a  normalized  value  at  the  center  of  the  plate  (Graph 
2).  Here,  the  first  clear  advantage  of  the  special  purpose  elements  manifests  it¬ 
self.  Element  QHl  only  approximates  the  traction-free  edges  condition  (20%  error) 
while  all  the  special  prupose  elements  exactly  satisfy  this  condition.  Among  the 
traction-free  elements,  PL20  follows  the  exact  solution  very  well,  while  PL21  tails 
off  at  the  special  element/regular  element  interface  and  the  curve  for  PL19  has  in¬ 
creasing  slope  rather  than  decreasing  slope  as  in  the  exact  solution. 

With  no  further  quantities  of  interest  in  Problem  1  and  too  small  a  basis  for 
making  any  broad  observations,  a  second  example  problem  will  be  considered. 

Shown  in  Figure  4  is  a  circular,  annular  plate  (inner  radius,  b;  outer  radius, a; 
thickness,  2h)  with  the  inner  edge  free  and  the  outer  edge  clamped.  It  is  subject 
to  a  transverse  load  of  magnitude  P  distributed  as  a  line  load  along  the  inner  edge 
such  that:  P*2irbQ^  C2lII.  Again,  by  sysnetry,  it  is  possible  to  model  a  quarter  of 
the  plate  in  an  FEM  analysis.  Only  one  mesh  size  is  necessary  (since  displacements 
were  shown  to  converge  in  Problem  1):  4  elements  in  the  radial  direction  and  6  ele- 
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Graph  t.  Problem  1 .  My  Along  CD  (Normalized) 
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Figure  4.  Problem  2.  Annular  Plate 


lo: 

ments  in  the  circumferential  direction.  The  traction-free  edge  elements  are  placed 

along  the  inner  edge  ACGE.  .Note,  this  problem  also  tests  the  performance  of  the 

elements  as  general  quadrilaterals  which  approximate  the  circular  edge  with  straight 

line  segments.  .All  results  are  presented  for  the  case:  a=12.0  in,  b=4.0  in,  0  * 

o 

1.0,  E=1.0  X  10^  psi,  h  =  0.5  in  and  v  =  0.3. 

Since  this  problem  is  axisymmetric,  a  closed  form  solution  can  be  readily  ob¬ 
tained,  and  various  quantities  of  interest  can  be  compared  with  numerical  solutions. 

In  Graph  3,  the  deflections  along  line  AB  are  plotted.  As  before,  all  the  elements 
predict  this  quantity  well.  Furthermore,  along  any  other  radial  line  (e.g.,  CD,  GH) , 
the  curves  should  not  change  and  this  is  verified  in  the  numerical  results.  The  ro¬ 
tation  about  the  y-axis,  0^,  along  line  AB  is  shown  in  Graph  4,  Notice,  that  traction- 
free  elements  PL19  and  PL20  underpredict  the  rotations  at  the  free  edge.  This  is  un¬ 
desirable  in  itself  and  may  affect  the  predictions  of  moments  and  stresses  subse¬ 
quently.  ^ 

Graph  5  shows  the  variation  of  along  line  AB.  starts  as  a  positive  quanti¬ 
ty,  passes  through  zero,  and  is  negative  at  the  clamped  edge.  Observe  that  PL19  and 
PL20  behave  poorly  in  terms  of  predicting  the  moment  at  the  free  edge  and  at  the 
special  element/regular  element  interface.  Element  PL21  does  well  at  the  free  edge 
but  shows  slight  difficulty  at  the  special  element/reguiar  element  interface.  On 
the  other  hand,  the  non-traction-free  element,  QHl,  follows  the  analytic  curve  very 
closely  at  all  points. 

The  other  moment  of  interest,  along  line  AB,  is  plotted  in  Graph  6.  This 
moment  is  zero  at  the  free  edge  and  increases  to  a  maximum  negative  value  at  the 
clamped  edge.  Here,  all  the  elements  predict  the  moment  well  at  the  free  edge,  but 
at  the  special  element/regular  element  interface  all  the  special  purpose  elements 
are  poor.  PL19  and  PL20  grossly  overestimate  the  moments  (over  300%  error)  while 
PL21  exhibits  an  error  of  about  50%  at  the  interface.  .Also,  PL19  and  PL20  cause 
the  regular  element  at  the  interface  to  predict  the  moment  poorly.  Finally,  at  the 
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Graph  3.  Problem  2.  Deflection  Along  AB 
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Graph  4  Problem  2.  Rotation  Along  AB 
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Graph  6.  Problem  2.  Along  AB 
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clamped  edge  (where  element  QHl  is  always  used) ,  the  moment,  M^,  is  underestimated 
in  all  cases.  More  will  be  said  on  this  shortly. 

In  assumed -displacoaent  elements,  much  work  has  been  done  to  determine  those  lo¬ 
cations  within  an  element  at  which  the  stress  (and  moment)  predictions  are  best 
[23,24,16].  These  points  are  called  "optimal  sampling"  points  and  depend  on  the  or¬ 
der  of  the  polynomials  used  in  the  displacement  interpolations  for  the  elements. 
Though  such  points  have  not  been  rigorously  defined  for  hybrid-stress  elements,  it 
is  at  least  reasonable  to  assume  that  these  points  will  not  lie  on  the  edges  of 
these  elements.  Therefore,  in  Problem  2,  line  AB  (Figure  4)  is  not  expected  to  be 
a  line  of  "optimal  sampling".  Moreover,  if  the  "optimal  sampling"  behavior  of  hy¬ 
brid-stress  elements  is  like  that  in  assumed -displacement  elements,  then  for  8 -node 
elements,  the  "optimal  sampling"  points  are  the  points  used  for  2x2  Gauss  Quadra¬ 
ture. 

Referring  to  Figure  4,  line  CD  is  a  typical  radial  line  which  contains  the  2x2 
Gauss  points.  In  Graphs  7  and  8,  the  normal  and  targential  moments,  respectively, 
along  line  CD  are  plotted.  Notice  that  even  along  a  line  of  "optimal  sampling"  trac¬ 
tion-free  elements  PL19  and  PL20  are  unable  to  predict  moments  well.  The  large  mo¬ 
ment  at  the  free  edge  in  Graph  7  is  largely  underestimated  by  PL19  and  PL20,  while 
at  the  special  element/regular  element  interface  both  moments  are  predicted  poorly. 
On  the  other  hand,  special  purpose  element  PL21  2md  regular  element  QHl  both  follow 
the  analytic  solutions  remarkably  well  along  this  line;  though  QHl  still  only  approx 
imates  the  traction-free  condition  since  the  tangential  moment  is  not  exactly  zero 
at  the  free  edge  (Graph  8). 

Recall,  that  in  Graph  6,  the  tangential  moment  was  underestimated  at  the  clamped 
edge.  Along  a  line  of  "optimal  sampling"  this  discrepancy  vanishes;  as  shown  in 
Graph  8  this  moment  is  predicted  adequately  at  the  clamped  edge  in  all  cases.  This 
gives  further  credibility  to  the  idea  of  evaluating  stresses  at  "optimal  sampling" 
points  rather  than  arbitrary  or  convenient  points  in  an  element. 
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Graph  7.  Problem  2.  Normal  Moment  (Mp)  Along  CD 
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Graph  8.  Problem  2.  Tangential  Moment  (Mt)  Along  CD 
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Based  on  these  observations,  element  PL21  seems  the  best  candidate  to  use  as  a 
single-layer,  traction-free  edge  plate  element.  Unlike  PL19  and  PL20,  it  predicts 
both  rotations  and  displacements  well.  Though  it  predicts  moments  better  than  PL19 
and  PL20  in  general;  along  a  line  of  “optimal  sampling"  it  predicts  moments  ex¬ 
tremely  well  and  shows  a  clear  advantage  over  those  elements.  As  far  as  element 
QHl  is  concerned,  ignoring  its  inherent  deficiency  of  approximating  the  traction- 
free  edge  condition,  it  performs  extremely  well  in  all  cases. 


Ill 


4 .  Summary  and  Concluding  Remarks 

In  this  study,  a  single-layer,  isotropic  plate  element  with  a  straight  traction- 
free  edge  has  been  developed  based  on  the  hybrid-stress  model.  It  can  be  coupled 
with  the  regular  version  of  this  element  (element  QHl  Cl33)  to  perform  the  numer¬ 
ical  analysis  of  isotropic,  arbitrary  thin  to  moderately  thick  plates  with  trac¬ 
tion-free  edges. 

The  element  assumes  a  Mindlin-type  displacement  behavior  which  not  only  results 
in  displacement  interpolations  that  are  only  required  to  be  C°  continuous  but  fa¬ 
cilitates  satisfying  the  3-D  equilibrium  equations  and  the  free-surface  conditions 
of  the  element.  By  including  all  components  of  stress,  the  transverse  shear  stresses 
are  present  for  the  analysis  of  moderately  thick  plates  where  their  contributions 
are  significant.  The  stress  assximption  for  the  element  is  chosen  so  that  along  one 
of  its  edges  the  tractions  (i.e.  appropriate  stresses)  are  zero;  and,  the  B-q  re¬ 
lation,  a  necessary  condition  in  the  assumed-stress  hybrid  formulation,  is  satis¬ 
fied.  The  element  mapping  is  such  that  it  can  take  on  a  general,  quadrilateral 
shape. 

Though  the  number  of  stress  fields  which  satisfy  these  conditions  is  theoreti¬ 
cally  limitless,  arguements  of  minimizing  the  number  of  stress  parameters  quickly 
narrows  the  choices.  Extensive  numerical  experimentation  and  testing  results  in 
the  conclusion  that,  of  the  numerous  stress  fields  considered,  the  21-6  stress 
field  of  element  PL21  is  the  best  for  the  applications  considered. 

Element  PL21  yields  a  stiffness  matrix  which  exhibits  correct  rank,  avoids  "lock¬ 
ing"  in  the  thin  plate  limit,  and  yields  displacements  which  converge  as  the  num¬ 
ber  of  degrees  of  freedom  increases.  Moreover,  in  a  numerical  analysis,  element 
PL21  predicts  displacements,  rotations,  and  stresses  (moments)  quite  well  in  gen¬ 
eral,  and  especially  well  along  lines  containing  the  2x2  Gauss  points  which  appear 
to  be  "optimal  sampling"  points  from  which  to  obtainor  extrapolate  stress  (moment) 
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information. 

However,  it  is  important  to  realize  that  in  all  cases  element  QHl  performs  well 
in  predicting  displacements  and  stresses;  and,  its  only  shortcoming  is  its  inher¬ 
ent  inability  to  satisfy  the  traction-free  condition  exactly  along  an  edge.  This 
seems  to  indicate  that  in  single-layer,  isotropic  plates,  the  traction-free  edge 
condition  is  not  crucial  when  trying  to  obtain  a  numerical  solution  and  that  the 
extra  effort  to  develop  and'xise  special  elements  to  better  model  this  condition 
may  not  be  necessary. 

On  the  other  hand,  the  stresses  which  dominate  the  solution  of  multi-layer,  lam¬ 
inated,  composite  plates  (i.e.  have  minimal  effects  in  single-layer, 

isotropic  plates.  In  fact,  these  stresses  have  not  even  been  considered  in  any 
great  length  in  this  study.  For  this  reason,  the  advantages  of  the  traction-fcee 
edge  elements  over  element  QHl  may  not  be  evident  in  single-layer  plates. 

Therefore,  a  multi-layer,  traction-free  edge,  plate  element  based  on  elonent 
PL21  should  be  developed  and  tested.  This  element  when  compared  to  the  multi-layer 
version  of  QHl  should  prove  to  perform  better  and  exhibit  obvious  advantages  when 
analyzing  laminated,  composite  plate  structures  with  traction-free  edges. 
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CHAPTER  4 

ELASTIC-PLASTIC  ANALYSIS  OF 
SINGLE -LAYER  PLATES 

ABSTRACT 


Two  alternate  hybrid-stress -based  functionals  are  examined  for  the  incremental 
elastic-plastic  static  analysis  of  single  layer  plates.  Material  nonlinear  effects 
are  incorporated  via  the  initial -stress  approach  so  that  an  equivalent  nodal  force 
vector  is  defined  and  the  stiffness  remains  constant  throughout  the  incremental 
loading.  The  alternate  functionals  differ  in  the  incremental  stress  which  is  assumed 
to  satisfy  equilibrium;  in  the  first,  it  is  the  actual  stress  increment,  and  in  the 
second  it  is  the  elastic  stress  increment.  Results  are  presented  for  two  example 
problems,  and  comparisons  of  the  alternate  functionals  and  plausible  iteration 
schemes  are  given.  The  effects  of  variation  of  pertinent  solution  parameters  are 
also  shown.  A  4-node  hybrid-stress  plate  element  based  on  a  Mindlin-type  displace¬ 
ment  field  is  used  for  most  cases;  however,  limited  results  are  also  presented 
using  an  8-node  plate  element,  thus  permitting  comparisons  of  the  relative  effi- 
ciences  of  the  tswo  elements. 
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1.  Introduction 

In  the  elastic-plastic  analysis  of  structures  by  the  finite-element  method,  three 
incremental  procedures  are  most  commonly  used;  the  tangent  stiffness,  initial-strain, 
and  initial-stress  methods.  In  the  first,  nonlinear  effects  axe  incorporated  via  an 
updated  stiffness,  whereas  in  the  latter  two  methods,  the  structxire  stiffness  re¬ 
mains  constant  and  plasticity  effects  are  included  as  equivalent  nodal  loads.  In 
most  cases,  the  analyses  are  based  on  assumed-displacement  formulations;  textbook 
accounts  and  appropriate  references  may  be  found,  for  example,  in  References  1  and 

■7 

A  viable  alternative  to  the  assumed-displacement  model  for  both  linear  and  non¬ 
linear  analyses  is  the  hybrid-stress  nodel.  This  model  is  based  on  a  modified  com¬ 
plementary  energy  principle.  Compatible  boundary  displacements  and  equilibrating 
intraelement  stresses  are  independently  interpolated;  the  stress  parameters  are 
eliminated  on  the  element  level  resulting  in  a  conventional  element  stiffness  matrix. 
Elastic-plastic  analyses  based  on  the  hybrid-stress  model  D-lOH  have  been  reported  by 
Yamanda  et  al  Csl,  Luk  C^],  Spilker  and  Pian  C5,7!],  Horrigmoe  and  Eidsheim  [l83, 
and  Barnard  and  Sharman  C9!];  these  include  applications  of  all  three  elastic-plaistic 
procedures.  A  sxirvey  of  incremental  hybrid-stress  formulations  for  nonlinear  prob¬ 
lems  has  been  presented  by  Pian  Cl03. 

The  initial-stress  procedure  was  first  used  in  conjunction  with  the  assumed-dis¬ 
placement  model  by  Zienkiewicz  et  al  Cll3-  In  this  approach,  the  effects  of  mater- 

f 

ial  nonlinearity  aire  accounted  for  by  a  fictitious  initial  stress,  a  ,  equal  to  the 

difference  between  an  assumed  elastic  stress  increment,  hc* ,  and  the  actual  stress 

increment,  Aa  (see  Figure  1).  The  equivalent  load  vector  is  calculated  by  a  weighted 
-  ep 

f 

integral  of  a  ,  and  therefore  the  accuracy  of  the  numerical  scheme  will  be  dependent 
on  the  accuracy  of  intraelement  stress  distributions.  Since  the  hybrid-stress  model 
yields  improved  intraelement  stress  predictions  compared  with  analagous  assumed- 
displacement  elements  in  many  cases  [e.g.  12,13j,  this  model  appears  to  be  ideally 
suited  for  use  in  conjunction  with  the  initial-stress  approach. 


UNIAXIAL  STRAIN.  £ 


Figure  1 


Schematic  representation  and  definition  of  stress  and  strain 
quantities  for  the  initial-stress  approach  for  a  1-D  problem, 
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Alternate  hybrid-stress  functionals  using  the  initial-stress  approach  have  been 

derived  in  References  5  and  7.  These  functionals,  denoted  bv  and  II ,  differ 

'me  me  ’ 

in  the  stress  increment  assumed  to  satisfy  equilibrium;  in  it  is  the  (approx¬ 
imate)  actual  stress  increment  (Aa  in  Figure  1) ,  while  in  it  is  the  elastic 

stress  increment,  Ac"".  Extensive  numerical  studies  for  axisymmetric  structures, 
also  including  equilibrium  imbalance  corrections,  suggest  that  II is  the  more 
.  accurate  and  efficient  approach  Cs, 

However,  a  potential  disadvaintage  of  is  the  need  for  intraelemant  compatible 

displacement  interpolations  (not  required  in  II^^) .  For  plate  elements  based  on 

classical  thin  plate  theory,  the  required  continuity  intraelement  displacement 

fields  are  not  easily  constructed  and  a  II^^  approach  would  therefore  be  intractable 

Applications  of  the  hybrid-stress  model  to  elastic-plastic  plate  bending,  using 

classical  plate  theory  elements,  have  been  reported  in  References  8  and  9.  In  these 
I 

approaches,  as  in  ,  displacements  are  interpolated  only  on  the  element  boundar¬ 
ies  and  thus  no  difficulties  are  encountered  in  defining  the  continuity  inter¬ 
polations. 

Recently,  a  family  of  hybrid-stress  plate  elements  have  been  developed  for  which 
independent  transverse  displacement  and  rotations  are  assumed,  and  in  which  all  com 
ponents  of  stress  are  included  L14-16J;  the  elements  are  thus  applicable  for  thin 
and  moderately  thick  plates.  Because  only  C°  displacement  continuity  is  required, 
intraeleraent  displaceraents/rotation  interpolations  are  easily  constructed,  and 
thus  the  approach  is  possible.  The  family  of  elements  are  based  on  4-node, 
8-node,  and  12-node  Serendipity  shape  functions.  In  each  case,  the  element  stiff¬ 
ness  is  of  correct  rank,  and  the  plate  thickness  may  be  taken  arbitrarily  small 
without  inducing  solution  'locking'.  In  numerical  comparisons,  these  elements,  in 
general,  yield  superior  displacement  and  intraelement  stresses  in  comparison  with 
analagous  reduced  or  selectively  reduced  integration  assumed-displacement  Mindlin 
plate  elements.  Because  the  hybrid-stress  plate  elements  do  not  require  a  formal 
separation  of  flexural  and  transverse  shear  stiffness  contributions,  the 
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application  of  these  elements  to  nonlinear  problems  will  be  straightforward  and 
no  special  modifications  will  be  required  in  order  to  analyze  plates  from  moderate 
to  arbitrarily  small  thickness. 

In  the  present  study,  the  plate  elements  of  Cl'ijIS]  are  used  with  and 

me  me 

for  elastic-plastic  analysis.  Equilibrium  imbalance  corrections  are  included  in 
all  cases,  and  the  effects  of  solution  refinement  and  integration  sampling  points 
are  explored,  primarily  using  and  a  4-node  plate  element.  Comparisons  of 
and  [using  alternate  iteration  schemes)  are  presented  which  again  show  to 
be  the  better  approach.  Also,  limited  comparisons  of  the  4-node  and  8-node  elements 
are  presented  to  assess  the  relative  efficiency  and  accuracy  of  these  elements  in 
nonlinear  analysis. 


2.  Hybrid-Stress  Functionals  and  Matrix  Formulation 
Details  of  the  derivation  of  the  alternate  hybrid  functionals  from  the  incremen¬ 
tal  virtual  work  expression  are  found  in  Reference  5;  only  the  final  forms  and  ele¬ 
ment  matrix  definitions  will  be  given  here  for  the  sake  of  completeness.  In  vector 

form,  and  II may  be  expressed  as: 
me  me 


nl  Cia.Au) 
me  --  '• 


(Aa  +  a^)  S  (Ac  +  dV  - 
Au  dS  * 


A  £ 


dV 


CD 


S 


a 


n 


ca5:au) 


where 


+ 


Acr*dV 


* 

n 


A  e  dV 


(2) 


■  volume  of  the  nth  element, 
n 

S(7  *  portion  of  the  boundary  of  the  nth  element  over  which  tractions 

n 

are  prescribed. 
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S  »  material  elastic  compliance  matrix  (  e  =Sa) , 

c'’ ■  elastic  stress  vector 
f 

a  3  fictitious  initial  stress  vector  Ca  prescribed  quantity) . 

f 

a  =  actual  stress  vector  (i.e.  Ac  =  Aa  -  ) . 

u  -  displacement  vector 

£  3  "strains"  calculated  from  displacements,  u,  via  the  linear  strain- 
displacement  relations. 

T  3  prescribed  tractions. 

A(  )  3  increment  in  the  quantity  (  ). 


C  )' 


total  value  of  the  quantity  C  ) • 


Tne  term  R*  appearing  in  both  functionals  corresponds  to  the  equilibrium  imbalance 
correction  [  S]  and  is  given  by  (for  no  body  forces): 


a..-/ 


^oT  .  A 
a  A  e  dV 


-  /  t”’ 


Au  dS 


'^n  C3) 

If  the  total  stress  at  the  beginning  of  an  increment  satisfies  equilibrium,  and 
the  total  tractions  satisfy  the  mechanical  boundary  conditions,  R*  vanishes. 

In  ^^^1  the  actual  stress  increment,  Aa,  must  satisfy  the  equilibriian  equation; 


E  Aa  3  E(Aa'  -  a  ) 


(4) 


where  E  is  the  differential  operator  matrix  corresponding  to  the  homogeneous  equili¬ 
brium  equations.  In  ,  the  elastic  stress  increment,  Aa' ,  must  satisfy  equilibrium 


E  Aa  3  E(Aa^p  +  a^) 


(S) 


where  is  the  actual  elastic-plastic  stress  increment  (see  Figure  1).  In  view 
of  equations  (4)  and  (5) ,  appears  to  be  the  more  "consistent"  hybrid-stress 
approach,  and  may  be  viewed  as  a  modified  Hellinger-Reissner  approach  since  only 
a  portion  of  the  actual  stress  increment  satisfies  equilibrium.  It  should  also  be 
noted  that  the  second  integral  in  appears  as  a  surface  integral  over  the  ele- 
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aent  boundary  in  the  original  formulation  l5,73  (by  applying  the  Divergence  theor- 
om  to  this  term,  with  equation  so  that  displacements  need  be  defined  only  on 

the  element  boundary  (if.  is  elminated] .  However,  for  the  present  application 
where  intraelement  displacements  are  easily  interpolated,  the  present  form  is  pre¬ 
ferred.  N'ote  also  that  an  intraelement  displacement  field  is  always  required  when 

using  ,  and  when  R*  is  retained  in  . 
me  ’  n  me 

In  the  element  formal at i'ons,  stresses  are  expressed  in  terms  of  stress  parameters, 
3  (usually  in  polynomial  form); 

Aa  =  P  AS 
,  Aa »  P  AS 

such  that  the  appropriate  equilibrium  equations  (equations  (4)  and  (S))  are  exactly 
satisfied.  The  intraelement  displacements,  u,  are  interpolated  in  terms  of  nodal 
displacements,  q,  such  that  appropriate  interelement  continuity  is  guaranteed; 


for  H' 


for  H 


me 

II 

me 


C&a) 


(6b) 


Au  =  N  Aq 

The  linear  strain-displacement  relations  are  applied  to  equation  (7)  to  give 


A'^  »  B  Aq 


Equations  (6)  through  (8)  are  substituted  into  equations  (1)  and  (2),  and  the 
following  element  matrices  are  defined 


(8) 


H./ 


P^SP  dV 


(9a) 


/  p’'b 


G  ■  /  P‘B  dV 
V 


(9b) 


•  f  P^So 


Sa^  dV 


(9c) 
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<T  dV 


(Qd) 


>  f  .T-.o 


a  dV 
•>» 


(9e) 


'iQ 


/ 

•4  Af 


AT  dS 


C9f) 


/  /j 


T°  dS 


C9g) 


Following  CS,73  ,  the  stress  parameters  are  eliminated  on  the  element  level  from 
the  resulting  functionals  according  to; 


AS  *  H“^G  Aq  -  H‘^F^  for 

me 


(10a) 


AS  *>  H’^G  Aq 


for  H 


II 

me 


ClOb) 


When  equations  (10)  are  substituted  back  into  the  functionals,  both  and  II may 
be  put  in  the  form; 


.1,11 


me 


Z  [^1/2  Aq\  Aq  -  Aq’''(AQ  +  +  R°)] 


(11) 


where 

T  -1 

k  »  G  H  G  a  element  elastic  stiffness  matrix.  (12a) 

AQ  a  incremental  element  external  load  vector.  (12b) 

R°  a  Q°-F°  a  equivalent  element  load  vector  corresponding  to 

the  equilibrium  imbalance  correction.  (12c) 

f 

The  vector  Q  is  the  equivalent  element  load  vector  corresponding  to  the  fictitious 
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initial  stress  and  differs  for  the  two  functionals; 


*  g'^h'^F^  for 
2^  ~  ~  ~  me 


Q  a  F 


for  n 


II 

me 


C13a) 

(13b) 


Following  the  usual  assembly  operations,  the  stationary  condition  of  yields 
the  following  system  of  equations,  written. for  the  ith  increment  (fron  applied  load 
i-1  to  applied  load  i) : 


(14) 


where  the  starred  (  )•  quantitites  refer  to  the  assembed  system  and  K  is  the  global 
elastic  stiffness  matrix. 

To  solve  jsquations  (14),  K  is  factored  prior  to  the  first  load  increment,  after 

* 

which  Aq^  can  be  calculated  at  each  increment  by  the  forward/backward  substitution 

f* 

operations.  Note  that  is  not  known  for  the  ith  increment,  but  can  be  extrapo¬ 
lated  or  estimated  from  data  at  the  previous  increment.  In  the  present  applications 
f*  f* 

is  used  in  conjunction  with  iteration  within  a  load  increment.  The  equi- 
o  * 

libriimi  correction  vector,  ,  is  calculated  from  the  total  stresses  and  total 
external  load  at  the  beginning  of  an  increment,  which  are  known. 


3.  Elastic-Plastic  Material  Relations 

The  effects  of  material  nonlinearity  are  incorporated  in  the  fictitious  initial 
f  f 

stress,  a  .  In  order  to  calculate  o  ,  the  correct  "elastic-plastic"  stress  incre¬ 
ment,  Aa  ,  corresponding  to  an  increment  in  total  strain,  A  e  ,  is  required  (Figure 
-  ep  -. 

1) .  This  relation  is  given  by 


D  A  e 
-ep  - 


(15) 


where  is 
-ep 


termed  the  elastic-plastic  material  matrix. 


-Jk. 

i 
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The  matrix  is  determined  from  the  flow  theory  of  plasticity,  assuming  a 
yield  surface  FCa,<)  »  0  and  an  associated  flow  rule.  It  can  be  shown  that  [e.g.  3, 


iBiy  p  [  *(ifJ  p  (I)  ] 


where  D  is  the  elasticity  matrix  (a  »D  e )  and  is  the  slope  of  the  equivalent 
stress  versus  equivalent  plastic  strain  curve.  is  symmetric  and  positive  def¬ 
inite  and  equation  (16)  is  valid  also  for  elastic,  perfectly-plastic  materials 

p 

(i.e.  =0).  Note  also  that  D  is  stress-history  dependent  and  must  therefore 

t  ~  ep  • 

be  evaluated  for  each  sampling  station  at  each  incremental  step. 

In  the  present  analysis,  the  Huber -Mises-Hencky  initial  yield  criterion  is  adopted. 
For  3-D  stress  states  (as  in  the  plate  elements  to  be  used  here),  the  yield  surface 


F(a,K)  = 


“  -  a  )^+  (a  -  a  )^  +  (a  -  a  )' 

L  x  y  X  2  ^  y  z-’ 

2  2  2  1  — 

+  6(7  +6a  +60  I  -(7a( 

xy  xc  y2  J  0 


and  therefore 


3F  _ 

^0  C! 


1  -172  -1/2 
1/2  1  -1/2 
1/2  -1/2  1 


where  0  is  the  uniaxial  yield  stress. 

0 

The  examples  considered  herein  assume  elastic,  perfectly-plastic  material  behavior. 
For  strain  hardening  materials,  various  mathematical  models  Ce.g.  17-19j  could  be 
incorporated . 


4.  Calculation  of  Equivalent  Element  Loads 


The  expressions  for  the  equivalent  element  loads,  and  F°,  require  the  dis- 
^  0 

tributions  of  r  and  J  which,  in  the  plastic  range,  can  onl/  be  determined  nua- 
erically.  Thus,  and  are  evaluated  using  a  numerical  integration  rule  (Gauss 
quadrature  in  the  present  case)  so  that  and  (and  appropriate  matrices)  are 
evaluated  at  numerical  integration  stations.  A  flow  chart  indicating  the  steps' re¬ 
quired  for  and  F®  for  an  element  and  for  both  II and  is  given  in  Figure  2. 
Note  that  steps  unique  to  or  have  been  prefaced  by  "I”  or  "II",  respectively; 
otherwise  the  steps  are  identical  for  both  functionals. 

A  comtarison  of  the  two  precedures  shows  that  recuires  more  core  storage 

me 

(both  H"^G  and  must  be  stored  for  each  element)  and  more  operations.  As  a  rough 

benchmark  estimate  of  the  relative  computation  times  required  for  the  evaluation 

of  and  F°  for  and  3^^  (for  one  element),  the  number  of  multiplications 

required  for  the  matrLx  operations  in  Figure  2  can  be  determined.  Table  1  gives  the 

multiplication  counts  for  each  of  the  operations  in  Figure  2;  note  that  it  has  been 

assumed  that  all  matrix  operations  are  done  in  full  without  accounting  for  symmetry 

or  zeroes.  Summing  the  contributions,  and  ignoring  the  equilibriian  correction,  the 

total  multiolications  (for  one  element),  and  for  3^  amd  3^^  ,  respectively, 

me  me 

are: 


M  a  jn-n  +  n  n  n  n 
3  q  X  y  z  s 


M  «  n.n  +  n  n  n  n 
3  q  X  y  z  s 


I  "s  "3  *  *  2n^  *  Hg)!. 

|n^  >  ng  >  a(3n^  ^  2n^  >  n^)  | 


(19a) 


(19b) 


where  a  is  the  ratio  of  the  number  of  integration  stations  per  element  at  which 
yielding  occurs  to  the  total  nximber  of  integration  stations  per  element.  The 
remaining  parameters  have  been  defined  in  Table  1.  These  expressions  will  be  used 
later  to  compare  various  solution  strategies.  Further  discussion  of  the  alternate 
functionals  and  solution  procedures  is  found  in  Reference  [5] . 


26 


Figure  2.  Flow  chart  for  the  calculation  of  equivalent  loads  for 
an  element.  Note  that  subscript  t  implies  a  temporary 
quantity. 


At  each  integration 
station 


Table  1.  Benchmark  multiplication  counts  per  element  for  the 
evaluation  of  the  equi\‘alent  load  vector  correspon¬ 
ding  to  material  nonlinearity. 
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(1)  These  computations  are  necessary  only  if  yielding  occurs  at 
that  station. 

Legend : 

ttg  =  number  of  3's  per  element. 

n^  =  number  of  degrees  of  freedom  per  element. 

n^  a  number  of  stress  components. 

n^  =  number  of  z  integration  stations. 

n^  a  number  of  x  (or  integration  stations. 

n^,  ■  ntimber  of  y  (or  n)  integration  stations. 
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Finally  it  should  be  noted  that  the  operations  and  computation  time  when  using 
the  assumed-displacement  model  with  the  initial-stress  approach  Ullj  is  roughly 
equivalent  to  that  required  for  II CSj. 


S.  Iteration  Schemes 

Typically,  nonlinear  schemes  make  use  of  a  combination  of  loaui  increments  with 

iteration  within  each  load  increment  cycle.  In  the  present  study,  two  alternate 

★ 

iteration  schemes  are  used.  In  both,  equation  C14)  is  first  solved  for  /Iq^  cor¬ 
responding  to  an  applied  increment  in  external  load  and  including  any  fictitious 
forces,  ,  remaining  from  the  previous  step  and,  if  desired,  the  load  correspon¬ 
ding  to  equilibrixjm  imbalance.  In  the  first  iteration  scheme,  a  series  of  itera¬ 
tions,  governed  by  the  equation  (for  the  kth  iteration  within  the  ith  external 
loading  increment) . 


(20) 


are  performed  until  the  equivalent  load  corresponding  to  the  fictitious  initial 
stress  is  sufficiently  small;  ie.  until 


WT 


<  RCONV 


where  RCOMV  is  a  small  parameter  (e.g,  0.01),  ',(  )1  denotes  the  magnitude  (squared) 
of  the  vector  (  ) ,  and  corresponds  to  the  equivalent  load  vector  prior  to  the 
first  iteration  (competed  for  the  applied  external  load  increment).  This  scheme 
is  termed  iteration  scheme  A,  and  is  depicted.,  for  a  one-dimensional  stress-strain 
curve  in  Figure  3a.  Note  that  displacments,  stresses,  and  strains  are  continually 
updated  during  the  iteration  cycle.  IVhen  equation  (21)  is  satisfied,  the  solution 
proceeds  with  equation  (14)  for  the  next  increment  in  external  loading. 

Scheme  A  has  been  applied  to  both  II^  and  in  Reference 


Note  that  equation  (36)  in  Reference  Cll]  is  incorrectly  stated  and 
should  correspond  to  the  present  equation  (21) . 
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However,  it  is  believed  that  scheme  A  is  best  suited  for  I  ,  as  verified  in 

me 

the  studies  of  C5,7Z.  In  those  studies,  the  use  of  and  scheme  A  yielded  poor 

II  f 

predictions  of  strain.  From  Figure  3a,  it  is  apparent  that,  with  J  must  van¬ 

ish  in*  order  for  the  asstimed  stress  increment,  Aa7  to  be  equal  to  the  actual  stress 

increment,  Aa  .  In  addition,  from  equation  (5),  =  0  is  required  if  Ac  is  to 

-  ep  -  -  ep 

satisfy  equilibrium.  In  contrast,  from  Figure  1  and  equation  (4)  for  it  is  observed 

that  =  0  is  not  required  in  order  for  the  interpolated  stress  Ac  to  be  equal  to 

the  actual  stress,  Ac  ,  and  therefore  Ac  ^  to  satisfy  equilibrium. 

- ep  ~  ep  ^ 

•An  alternate  iteration  scheme,  termed  scheme  B  C5,7Z,  would  therefore  appear  to 
be  more  appropriate  for  In  this  scheme,  the  equation  governing  the  kth  itera¬ 

tion  within  the  ith  external  loading  increment  is 

K  Aq*  =  AQt  ^  qf , 

-  -k  -k-l  (-22) 

f* 

where  (i.e.  k^l)  is  zero.  This  scheme  (see  Figure  3b  for  an  illustrative  1-D 
case)  seeks  to  satisfy  the  condition 

-k  ~k-l  r2.'51 


which,  if  satisfied,  implies  that 

Ac^  =  Ac 

-1  ..ep.  f24) 

for  the  ith  external  loading  increment.  The  iteration  process  is  determined  to  be 


converged  when 


-  Cl/ 

hl'l 


RCONV 


is  satisfied.  In  scheme  B,  total  displacements,  stresses,  and  strains  are  updated 
only  after  equation  (25)  is  satisfied,  using  the  incremental  quantities  computed 
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(a)  LOAD  ITERATION  SCHEME  A 


Figure  3.  Schematic  1-D  representation  of  iteration  schemes  A  and  B. 


L 
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for  the  final  iteration  step. 

•An  advantage  of  scheme  B  with  is  that,  in  view  of  eq^lation  (24),  total  stresses 

may  be  updated  using  Aa.  Since  Aa  satisfies  equilibrium,  the  total  stress,  a®,  will 
be  in  equilibrium  and  no  equilibrium  imbalance  correction  is  needed.  If  R®  can  be 
ignored,  no  interior  displacement  field  need  be  defined  when  using  .  This  is 
particularly  important  for  elements  requiring  continuity,  but  is  of  no  real  con¬ 
sequence  for  C®  continuity  elements  (as  in  the  present  application). 

6.  Description  of' Plate  Elements 

The  plate  elements  to  be  used  are  taken  from  a  recently-developed  family  (4,3,  and 
12  node)  of  hybrid-stress  plate  elements  Cl4-16^.  The  elements  utilize  independent 
interpolations  for  the  transverse  displacement,  w,  and  cross-section  rotations,  9^ 
and  9y,  so  that  any  of  the  C®  continuity  families  of  shape  functions  may  be  used 
(in  these  cases,  the  Serendipity  shape  functions).  In  general,  all  components  of 
stress  are  included,  allowing  for  the  analysis  of  moderately- thick  plates. 

Results  presented  in  References  14-16  show  these  elements  to  yield  comparable  or 
superior  results  when  compared  to  the  analogous  assumed-displacement  plate  element 
(also  based  on  a  Mindlin-type  displacement  behavior  and  utilizing  reduced  or  selec¬ 
tive-reduced  integration).  Each  of  the  hybrid-stress  elements  of  this  family  has  a 
stiffness  of  correct  rank  (i.e.  no  spurious  zero  energy  modes)  and  the  elements  will 
not  'lock'  for  arbitrarily  thin  plates,  independent  of  machine  precision  used. 

Comparisons  made  in  Reference  16  suggest  that  the  8-node  and  12-node  elements  are 
the  more  accurate  per  degree-of-freedom  in  the  assembled  structure.  However,  12-node 
elements  may  be  impractical  for  general  applications.  Also,  for  nonlinear  analyses, 
where  computation  time  is  strongly  dependent  on  element-level  operations,  the  4-node 
element  may  be  preferred. . In  the  present  study,  both  the  4-  and  8-node  elements  will 
be  used  (with  most  emphasis  on  the  4-node  element) . 

The  elements,  both  of  which  allow  for  isoparametric  planforms,  are  shown  in  Fig¬ 
ure  4.  Briefly,  the  4-node  element  (termed  LH4  Cl43  ),  contains  12  degrees  of  free- 
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Figure  4.  Summary  of  plate  elements  UI4  and  Qlil . 
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dom  and  is  based  on  a  113  stress  field  (bilinear  x,y  variation  for  inplane  stresses). 
Note  that  a,  =  0  in  LH4.  The  8-node  element  (termed  QHl  ClS]  ),  contains  24  degrees- 
of-freedom  and  is  based  on  a  23S  stress  field  (bicubic  x,y  variation  for  inplane 
stresses)  with  #  0. 

It  is  important  to  observe  that  a  criticism  CsH  of  for  plate  analyses  -  -  - 
the  requirement  of  an  interior  displacement  field  yielding  interelement  contin¬ 
uity  -  -  -  is  no  longer  valid.  The  interior  displacement  fields  for  the  present 

plate  elements  (C°  continuity)  are  easily  constructed,  and  does  not  present 

*  me 

any  formulation  or  computational  difficulties. 

7.  Example  Problems  and  Numerical  Results 

The  problem  of  a  simply- supported  thin  beam  subjected  to  a  uniformly  distributed 

transverse  load  of  magnitude,  P,  has  been  chosen  as  a  first  test.  An  analytic  solu- 

tion  for  this  problem  is  available  C203  under  the  assumption  of  elastic,  perfectly- 

plastic  material  behavior.  Here,  the  beam  length  is  10.0  in,  depth  is  1.0  in., 

thickness  is  0.1  in.  The  elastic  material  constants  axe  E  =  10^  psi,  and  \>  =  0.3, 

.  .  4 

with  a  uniaxial  yield  stress  of  10  psi.  First  yielding  will  occur  at  the  beam 
center  at  a  load  of  1.33  psi,  and  the  fully-plastic  load  is  Ppp  =  2.0  psi. 

For  initial  comparisons,  the  4-node  plate  element,  LH4  [14], is  used  and  half  of 
the  beam  span  is  modelled  by  a  mesh  of  NDX  by  NDY  equal-sized  elements  (see  Fig¬ 
ure  5) . 

For  a  given  functional,  effects  of  various  parameters  can  be  examined;  these  in¬ 
clude  iteration  scheme  (A  or  B),  number  of  equal-size  load  increments  in  the  load 
range  from  initial  to  full  plasticity,  NINC  (i.e.  the  load  corresponding  to  first 
yield,  Pal.33  psi,  is  applied  first,  after  which  increments  in  load  corresponding 
to  AP» (2. 0-1. 33) /NINC  are  applied.  Loading  ceases  at  the  theoretical  fully  plastic 
load),  convergence  ratio,  RCONV,  maximijm  number  of  iteration  permitted  per  load 
step,  .MAXIT,  mesh  size  (NDX,  NDY),  and  integration  stations  for  computing  the  equiv¬ 
alent  loads  (NGX,NGY,NGZ) .  Unless  otherwise  stated,  NDX  ■  4,  NDY  *  1  is  used; 
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Figure  5.  Geometry,  mesh,  and  boundary  conditions  for  the  example  problem 
of  a  simply-supported  beam  (modeled  by  plate  elements)  subject 
to  uniform  load. 
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remaining  parameters  are  defined  on  each  figure.  In  all  cases,  the  predicted  nor¬ 
malized  center  transverse  displacement  (normalized  by  the  center  deflection 
at  first  yield)  versus  the  applied  load  (normalized  by  the  fully-plastic  load, 

Ppp)  is  compared  with  the  analytic  result.  In  a  first  series  of  tests,  only 

II will  be  used  (determined  to  be  the  superior  functional  in  C5,73).  In  all  cases 

me 

the  equilibrium  imbalance  correction  is  included  and  iteration  scheme  A  is  used. 

The  effects  of  NINC  (ntimber  of  load  steps  from  first  to  full  plasticity)  is  shown 
in  Figure  6  where  MAXIT  *  1  so  that  no  iteration  is  used.  Decreasing  load  step 
sizes  lead  to  an  improved  response,  particularly  near  the  fully-plastic  load.  The 
effects  of  maximum  permitted  iterations  per  load  step,  MAXIT,  for  load  increments 
NINC  =  8,4,  and  2  are  shown  in  Figures  7a  through  7c,  respectively.  In  each  case, 
MAXIT  =  1  (no  iteration),  10  and  20  are  used.  For  smaller  load  steps  (NINC  *  8, 
Figure  7a),  iteration  has  a  major  effect  only  near  the  fully  plastic  load;  however, 
increasing  MAXIT  from  10  to  20  has  no  effect  on  the  solution.  For  increasing  load 
step  size  (NINC  »  4,  Figure  7b),  MAXIT  has  a  greater  effect  and  differences  are 
observed  between  MAXIT  =  10  and  20.  For  large  load  steps  (NINC  *  2,  Figure  7c), 
solutions  are  poor  even  when  MAXIT  »  20. 

The  effects  of  convergence  ratio,  RCONV,  are  shown  in  Figure  8  for  NINC  *  4 
and  MAXIT  =  20;  a  value  of  RCONV  =  .01  has  been  used  previously.  Increasing  RCONV 

to  .1  has  a  discemable  effect  at  all  load  levels,  whereas  decreasing  RCONV  to 

.001  produces  significant  effects  only  for  the  last  load  step.  Arbitrary 
decreases  in  RCONV  without  a  corresponding  increase  in  MAXIT  will  not,  in  general, 
lead  to  a  significantly  improved  solution  as  MaXIT  will  terminate  the  iteration 

cycle  for  load  steps  in  severely  nonlinear  regions.  For  practical  purposes,  values 

of  RCONV  a  .01  and  MAXIT  *  10  are  believed  to  be  adequate.  Converging  results  can 
then  be  sought  by  increasing  the  nvnnber  of  load  steps  (smaller  load  increments) . 

The  computation  time  required  to  compute  the  equivalent  load  vector  is  nearly 
proportional  to  the  number  of  integration  stations  employed  (see  Table  1  and  equa¬ 
tion  (19b)).  The  question  of  through- thickness  and  inplane  stations  can,  however. 


Effect  of  maximum  permitted  iterations  per  load  step  for 
the  beam  problem. 
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I'igure  7  (concluded) 
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be  examined  independently.  The  through-thickness  distributions  used  for  all  stress 

components  in  the  present  hybrid-stress  plate  elements  should  correspond  to  the 

'exact'  distributions  for  thin  to  moderately-thick  plates  in  the  linear  elastic 

region.  Yielding  will,  in  general,  initiate  at  and  progress  from  the  plate  upper/ 

lower  surfaces,  resulting  in  discontinuous  through-thickness  stress  distributions 

(for  a  ,  a  ,  and  a  ).  This,  in  principle,  requires  that  increasing  nvimbers  of  z 
X  y  xy 

integration  stations  be  used  (which  therefore  locates  stations  closer  to  the  upper/ 
lower  surfaces).  The  effects  of  NGZ  (number  of  z  stations)  are  shown  in  Figure  9a 
where  MGX=NGY=3  and  MINC=4.  Decreasing  NGZ  from  3  (used  in  all  previous  cases)  to 
2  severely  stiffens  the  solution,  whereas  increasing  from  NGZ=3  to  NGZ=4  produces 
only  a  slightly  improved  solution.  In  view  of  the  increased  computation  time  asso¬ 
ciated  with  NGZ=4  compared  with  NGZ=3  with  only  marginal  improvement  in  the  results 
obtained,  the  value  NGZ=3  would  appear  to  be  the  appropriate  choice. 

The  selection  of  the  mamber  of  inplane  stations,  .VGX  and  NGY,  is  dependent  on 
the  accuracy  of  the  predicted  x-y  intraelement  stress  distribution;  i.e.  there  is 
no  merit  in  sampling  stresses  (from  which  equivalent  loads  are  calculated)  at  lo¬ 
cations  where  predicted  stresses  may  be  severely  in  error.  For  assumed-displacement 
elements,  optimal  sampling  points  can  be  defined  Ce.g.21,22]  and  used.  No  such 
points  have  been  rigorously  defined  for  hybrid-stress  elements;  however,  stress 
results  given  in  Cl4]  for  linear  analyses  suggest  that  LH4  yields  good  intraele¬ 
ment  stress  predictions  throughout  the  element,  and  therefore  increasing  NGX  and 
NGY  should  produce  improved  results.  Figure  9b  shows  the  effects  of  NGX  and  NGY 
(for  NGZ  =  3) .  The  results  for  NGX=i'«jY=3  are  marginally  improved  compared  with 
NGX»NGY=2.  The  solution  for  NGX»NGY=1  is  slightly  stiffer  until  the  final  load 
step  where  an  excessively  flexible  solution  is  obtained;  once  this  point  is  fully 
yielded,  the  entire  element  is  assumed  to  be  yielded  and  a  mechanism  forms  pro¬ 
ducing  the  excessive  flexibility.  In  view  of  the  results  obtained  and  computation¬ 
al  considerations,  NGX»NGY=2  is  recommended. 

The  1  point  inplane  integration  can  lead  to  acceptable  results  if  a  more  refined 
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mesh  is  used,  as  shown  in  Figure  10.  Using  the  NDX=4,  NDY=1  case  as  a  basis  for 
comparison  (with  NGX=NGY=1) ,  the  excessive  flexibility  no  longer  appears  when  NDX* 

S,  NDY=1  is  used  (note  that  a  fully  plastic  zone  has  not  yet  been  reached  in  this 
solution  at  P»Ppp) .  Increasing  to  NDX=8,  NDY=2  shows  a  slight  improvement. 

Using  the  benchmark  multiplication  counts  given  in  Table  1,  total  multiplications 
(for  a  load  increment  or  iteration  cycle,  and  a  specified  number  of  elements,  as¬ 
suming  yielding  occurs  at  all  integration  stations)  have  been  computed  for  a  number 
of  cases  of  interest  and  are  given  in  Table  2.  Based  on  these  benchmark  counts, 
the  use  of  a  2  by  2  block  of  elements  with  NGX=NGY=1  in  each  element  should  re¬ 
quire  a  similar  computational  effort  compared  with  1  element  and  NGX=NGY=2.  This 
is  analogous  to  an  3x2  mesh  with  NGX*NGY=1  and  a  4x1  mesh  with  NGX*NGY®2.  Results 
for  these  two  cases,  shown  in  Figure  11,  indicate  that  the  first  approach  is  su¬ 
perior,  despite  the  fact  that  the  predicted  elastic  center  deflecion  is  nearly  iden¬ 
tical  for  both  meshes.  This  improvement  would  therefore  appear  to  correspond  to 
superior  stress  accuracy  at  the  element  centroid  compared  with  the  2x2  Gauss 
stations. 

The  results  presented  thus  far  have  been  obtained  using  .  Similar  studies  have 

me 

been  carried  out  using  which  verify  the  general  trends  (with  regard  to  NINC, 
MAXIT,  RCON,  and  integration  stations)  observed  for  .  A  more  important  con¬ 
sideration  is  the  effect  of  iteration  scheme.  Figure  12a  shows  results  obtained 
using  iteration  scheme  A  for  a  coarse  load/iteration  solution  (NINC=4,  MAXIT=1) 
versus  a  refined  solution  (NINC=8,  MAXITalO).  In  both  cases,  NDXs4,  .NDY=1,  RC0NV= 
.01,  NGXaNGY®2,  NGZ*3  are  used.  The  coarse  solution  is  stiff  (as  expected)  whereas 
e.xcessive  flexibility  is  observed  for  the  refined  solution;  the  predicted  center 
displacement  exceeds  the  exact  value  for  P/Ppp  0.75.  Such  behavior  is  judged 
unacceptable  and  it  would  appear  that  iteration  scheme  A  is  not  well  suited  to 

n  ^  .  Note  that  further  solution  refinement  (results  not  shown)  leads  to  increased 
me 

flexibility. 
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Table  2.  Multiplication  count  comparisons  for  various  cases  of  interest. 


Functional 

"s 

\'"y 

^2 

%lemsP^ 

Mults. 

I 

LH4 

6 

2 

3 

1 

11,484 

I 

LH4 

• 

6 

1 

3 

1 

3,168 

I 

LH4 

3 

2 

3 

1 

2,592 

I 

LH4 

3 

1 

3 

1 

945 

II 

LH4 

6 

2 

3 

1 

10,860 

II 

LH4 

6 

1 

3 

1 

2,814 

II 

LH4 

3 

2 

3 

1 

2,256 

II 

LH4 

3 

1 

3 

r 

663 

II 

LH4 

6 

2 

4 

1 

14,436 

II 

LH4 

6 

1 

3 

4 

12,672 

II 

QHl 

6 

2 

3 

1 

13,008 

II 

QHl 

6 

3 

3 

1 

28,578 

II 

QHl 

3 

2 

3 

1 

3,540 

Notes:  (1]  "I"  and  "II"  correspond  to  and  ,  respectively. 

(2)  For  element  LH4,  n  =  12  and  n^  =  11. 

q  p 

For  element  QHl,  n  »  24  and  n-  *  22. 

q  p 

(3)  Number  of  elements  for  which  the  equivalent  load  is  com¬ 
puted.  It  is  assumed  that  yielding  occurs  at  all  Integra 
tion  stations. 


Jjgure  11,  Comparison  of  alternate  mesh/integration  rule  strategies  for 
the  beam  problem. 
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Figure  12.  Beam  example  results  obtained  using  hybrid  functional  I  and 
iteration  schemes  A  and  B. 


150 


z 

LU 

s 

UJ 

o 


CL 

CO 

O 

O 

LU 

N 

< 

IS 

cc 

o 

z 


0) 

3 


J 


(b)  Results  using  Iteration  Scheme  B 
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3y  comparison,  results  obtained  by  using  with  iteration  scheme  B  for  the 

coarse  and  refined  solutions  axe  shown  in  Figure  12b.  For  the  coarse  solution, 

scheme  B  yields  a  more  flexible  solution  compared  with  scheme  A,  whereas  the 

refined  solution  by  scheme  B  is  stiffer  than  that  of  scheme  A.  The  scheme  B 

refined  solution  is  in  good  agreement  with  the  analytic  solution,  although 

slightly  more  flexible  than  the  analytic  solution,  and  is  therefore  the  preferred 

scheme  with  3^  .  Note  that  further  refinement  of  the  scheme  B  solution  yields 
me 

no  discemable  change. 

A  comparison  of  with  scheme  B  and  3^^  with  scheme  A  for  coarse  and  refined 
solutions  is  shown  in  Fig\xre  13.  Recall  that  for  all  cases,  load  has  been  applied 
only  up  to  the  theoretical  fully-plastic  load,  Ppp*  Thus,  for  example,  the  re¬ 
fined  3^^  run,  which  underestimates  the  displacements,  would  produce,  effec- 
mc 

tively,  infinite  displacements  by  application  of  an  additional  small  increment 

in  external  load  beyond  P  =  P  . 

FP 

Both  fxinctionals  (with  an  appropriate  iteration  scheme)  appear  to  lead  to  con¬ 
vergence  to  the  analytic  solution.  Such  observations  were  also  made  in  References 
5  and  7  for  axisymmetric  structures.  However,  3^^  requires  less  computation  time 
per  cycle;  see  References  C5.73  and  benchmark  multiplication  counts  in  Table  2. 
Therefore,  3^^  would  again  appear  to  be  the  preferred  functional. 

Before  presenting  results  for  a  more  general  plate  problem,  several  observations 
related  to  computation  time  should  be  made.  In  the  cases  considered  here,  ail 
components  of  stress  have  been  included  (i.e.  n^  =  6) .  For  moderately  thick 
plates,  all  components  should  be  retained,  whereas  for  thin  plates  o  ,  a  and 
Cj  are  negligible  and  could  be  ignored  in  the  computation  of  equivalent  loads. 

From  Table  1,  total  multiplications  are  strongly  influenced  by  the  number  of 
stress  components  retained,  n^,  and  for  selected  cases  computed  in  Table  2,  it 
is  apparent  that  a  reduction  from  n^  *  6  to  n^  ■  3  can  lead  to  substantial  reductions 
in  multiplications  (i.e.  by  a  factor  of  4-5  times).  Results  (not  shown)  for  the 
beam  problems  (thin  plate)  show  the  two  solutions  (n^  *  3,6)  to  be  essentially  iden- 


agains 
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tical . 


Results  presented  also  suggest  that  NGZ=3  is  adequate  for  through-thickness 

integration.  If  the  plate  is  thin  so  that  n^  =  3  can  be  used,  the  remaining 

stresses  ,J  ,a  )  are  zero  at  2=0,  corresponding  to  the  center  integration 
X  y  xy 

station  for  NGZ=3,  amd  this  point  can  be  skipped  in  the  integration  loop.  The 
result  is  a  reduction  of  multiplications  by  approximately  33%.  Finally,  it  is 
noted  that  in  the  elastic-plastic  plate  analysis  of  Horrigmoe  and  Eidsheim  Csl, 
the  yield  surface  has  been  expressed  in  terms  of  moments,  thereby  avoiding  through¬ 
thickness  numerical  integration.  There  is  no  difficulty  in  adopting  such  a  pro¬ 
cedure  with  the  present  element  and  functionals  by  operating  with  analytically 
integrated  stresses.  However  this  approach  has  not  been  pursued  because  of  the 
inherent  smearing  of  actual  partial  yielding  effects  through  the  plate  thickness. 

The  second  example  problem  is  a  simply  supported  square  plate  of  side  length 
10.0  in.,  thickness  0.1  in.  subjected  to  uniform  load  of  magnitude  P  (Figure  14). 
The  material  is  assumed  to  be  elastic  perfectly  plastic  with  E  «  10^  psi,  v  = 

0.3,  and  =  10^  psi.  Upper/lower  bounds  for  this  problem  correspond  to  loads  of 
6.63  and  6.23  psi,  per  Hodge  and  Belytschko  [23].  Symmetry  is  utilized  and  a  quar¬ 
ter  of  the  plate  is  modelled  by  a  uniform  NDX  by  NDY  mesh. 

Because  is  the  computationally  preferred  scheme,  only  with  iteration 

scheme  A  and  including  equilibrium  imbalance  correction  is  used  in  this  example. 

In  view  of  the  parametric  studies  for  the  beam  problem,  the  parameters  MAXIT=10, 
RCONVs.Ol,  NGX=NGY=2,  and  NGZ=3  are  used  unless  otherwise  stated.  A  load  corres¬ 
ponding  to  P  =  3.47  psi  is  applied  first,  after  which  the  range  P  =  3.47  to  6.63 


psi  is  divided  into  NINC  equal  increments  in  load.  Results  are  presented  in  terms 


of  the  normalized  center  transverse  displacement  , (where  D  is  the  flex¬ 

ural  stiffness,  and  M  is  the  plastic  moment  at  full  yield)  versus  the  normalized 

2-  P 

pt  ^ 

total  load  f  -  . 

P 

The  effects  of  iteration  (using  a  2  by  2  element  mesh  for  which  the  elastic 


tip  deflection  is  in  error  by  less  than  2%)  are  shown  in  Figure  15. 


Boundary  Conditions: 

Sides  AB  and  AD:  w=0p,aO 
Side  BC:  dy-0 
Side  DC:  *0 

Loading:  Uniform  of  magnitude  P  (Ib/in) 

Mesh:  NDX  by  NDY  mesh  of  elements  LH4  (shown) 
or  QHI. 


Figure  14.  Mesh,  geometry,  and  boundary  conditions  for  the  example 
problem  of  a  square  plate  under  uniform  load. 


Q 

H* 

O 

z 

D 

u. 


^  o 

II  ¥ 

t  t 

X  X 

<  < 

2  s 


o 

ca 

II 

K 

X 

< 


9  □  O  < 

cc 

CQ 

> 

I 


cm  ^ 

P  U  II 

o  >>- 

>^(nZZ- 

Z  11  .  II 

o9nx  X 
o  =  c5o  D 
QC  zzzz 


in 

CM 


o 

CM 

d 


ure  15.  Effects  of  maximum  permitted  iterations  per  load  step  for  the 
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Increasing  MAXIT  yields  a  more  flexible  solution.  The  MAXIT*10  solution  will 
produce  a  fully-plastic  load  slightly  above  the  upper  bound  value,  whereas  the 
MAXIT=20  solution  will  fall  within  the  upper/ lower  bounds. 

The  effects  of  mesh  size  and  inplane  integration  stations  are  shown  in  Fig¬ 
ures  16a  and  16b.  The  comparison  between  the  case  of  a  2  by  2  mesh  with  NGX= 

NGY=2,  and  the  case  of  a  4  by  4  mesh  with  NGX=NGY=1  (approximately  equivalent 
in  terms  of  computational  effort).  Figure  16a,  suggests  that  the  latter  approach 
is  superior;  although  the  predicted  fully-plastic  loads  will  not  differ  substantially, 
details  of  the  load-deflection  results  are  clearly  different  beyond  a  normalized  load 
of  p>18.0.  The  effects  of  mesh  CNDX,.\DY),  for  NGX=NGY=1  (Figure  16b)  show  little 
change  in  predicted  response  beyond  a  3  by  3  mesh. 

In  Reference  lll6],  comparisons  of  the  hybrid-stress  family  of  plate  elements 
suggest  that  the  8-node  element,  QHl,  is  in  general  more  accurate  than  the  4-node 
element,  LH4,  per  degree  of  freedom  in  the  assembled  mesh  (an  exception  is  the 
simply-supported  plate  under  uniform  load).  A  reasonable  comparison  of  these  two 
elements  for  elastic-plastic  analysis  corresponds  to  a  2  by  2  block  of  LH4  elements 
with  NGX=NGY=1  replaced  by  one  QHl  element  with  NGX=NGY=2.  Benchmark  multiplication 
estimates  in  Table  2  suggest  that  the  LH4  analysis  will  require  slightly  less 
computational  effort.  Results  obtained  using  a  4x4  mesh  of  LH4  elements  with 
MGX=NGY=1,  compared  with  a  2x2  mesh  of  QHl  elements  with  NGX=NGY=2  are  showTi  in 
Figure  17.  The  two  schemes  are  found  to  produce  nearly  identical  results.  Actual 
CPU  times  show  that  the  LH4  analysis  requires  approximately  80%  of  the  time  re¬ 
quired  for  the  QHl  analysis.  It  should  be  noted  that  for  the  elastic  analysis  of 
a  simply-supported  uniform  loaded  plate,  LH4  is  more  efficient  than  QHl.  However 
for  other  boundary  conditions/ loading,  QHl  is  more  accurate  per  degree  of  freedom 
Cl6j  and  advantages  of  QHl  over  LH4  in  the  elastic-plastic  analysis  may  be  found. 
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Figure  16.  Effects  of  mesh  and  inplane  integration  rule  for  the  i)latc  problem. 
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8.  Concluding  Remarks 

The  elastic-plastic  analysis  of  plates  has  been  performed  by  using  alternate  hybrid- 
stress  functionals  based  on  the  initial-stress  approach  for  material  nonlineari¬ 
ties.  Accurate  results  have  been  obtained  using  both  with  iteration  scheme  B. 

me 

and  with  iteration  scheme  A.  In  principle,  has  the  potential  difficulty 

of  requiring  an  intraelement  displacement  field.  However,  this  poses  no  problems 
for  the  plate  elements  used  here  for  which  such  interpolations  are  easily  constructed. 
On  the  basis  of  computational  effort,  II is  preferred  over 

Results  obtained  suggest  that  a  2  by  2  inplane  and  3  point  through -thickness 

integration  rules  are  adequate  for  evaluation  of  the  equivalent  nodal  loads  corres¬ 
ponding  to  plasticity  effects  for  the  4-node  element.  Although  no  such  examples  have 
been  shown,  the  elements  and  procedure  used  here  can  also  be  applied  to  moderately- 
thick  plates.  Because  the  4-node  plate  element  produces  good  intraelement  stress 
distributions,  the  use  of  a  2  x  2  inplane  Gauss  rule  leads  to  improved  results 
compared  with  a  1  point  rule,  for  fixed  mesh  size.  However,  the  1  point  rule,  coupled 
with  a  more  refined  mesh  may  be  more  effective  from  a  computation  time  versus  accuracy 
viewpoint.  Significant  computational  savings  can  also  be  realized  for  thin  plates 
by  using  only  the  inplane  stresses  to  compute  the  equivalent  loads  and  bypassing 
the  z=0  integration  station  of  the  3-point  through-thickness  rule. 

A  single  comparison  between  die  4-node  element  and  8-node  element  CQHl) 

shows  that  the  use  of  the  simple  LH4  element  is  more  efficient  computationally. 

However,  for  other  loadings  and  boundary  conditions,  this  observation  could  be 
reversed.  In  QHl,  a  is  included,  whereas  a  =  0  in  LH4.  If,  for  example,  the  pres- 
ent  elements/analysis  were  extended  to  the  material  nonlinear  analysis  of  multi¬ 
layer  laminated  composite  plates,  a  QHl-type  element  should  be  preferred  because 
of  the  dominant  role  of  near  free-surfaces  for  relatively  thin  laminated  plates 
C24,  25]. 
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EDGE  SINGULARITIES  IN  .ANISOTROPIC  COMPOSITES 

^  ■■ 

ABSTRACT 

The  stress  singularity  at  the  vertex  of  an  anistropic  wedge  has  the 
fom  r''’F.'r,5)  as  r  — 0  where  0<^<1  and  F  is  a  real  function  of 
the  polar  coordinates  (r,5).  In  many  cases,  F  is  independent  of  r. 

The  explicit  form  of  F(r,9)  depends  on  the  eigenvalues  of  the  elasticity 
constants,  called  p  here,  and  on  the  order  of  singularity  ',\?hen  x 

is  real,  ?  =  <•  If  <  is  complex,  5  is  the  real  part  of  <.  The  p's 

are  all  complex  and  consist  of  5  pairs  of  complex  conjugates  which  reduce 
to  ±  i  when  the  material  is  isotropic.  The  function  F  depends  not  only 
on  p  and  <,  it  also  depends  on  whether  p  and  <  are  distinct  roots 
of  the  corresponding  determinants.  In  this  paper  we  present  the  function 

Ffr,9}  in  terms  of  p  and  <  for  the  cases  when  p  and  <  are  single 

roots  as  well  as  when  they  are  multiple  roots.  The  reiationsnip  between 

'  the  complex  variable  Z  introduced  in  the  analysis  and  the  oolar 

i 

1 

I  coordinates  (r,9)  is  interpreted  geometrically.  .After  presenting  the  form 

of  F  for  individual  cases,  a  general  form  of  F  is  given  in  Eq.  (64). 

'.Ve  also  show  that  the  stress  singularity  at  the  crack  tip  of  general 
anisotropic  materials  has  the  order  of  singularity  Z-h  which  is  a 
multiple  root  of  order  3.  The  implication  of  this  on  the  form  F(r,9)  is 


A 

♦ 

"4 


discussed. 
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1.  INTRODUCTION 

For  isotropic  materials,  use  of  the  biharmonic  function,  or  the  Airy 
stress  function,  seems  to  bfe  the  universal  approach  in  the  analysis  of 
stress  singularities  [1-4].  There  appears  to  be  no  universal  approach  in 
analyzing  the  stress  singularities  in  anisotropic  materials.  Lekhnitskii 
[5]  introduced  two  stress  functions  to  analyze  general  anisotropic  materials. 
His  approach  was  used  by  Wang  and  Choi  [6]  to  study  the  thermal  stresses 
at  the  interface  in  a  layered  composite.  Green  and  Zerna  [7]  employed  a 
complex  function  representation  of  the  solution.  Their  approach  was  used 
by  Bogy  [8]  and  Kuo  and  Bogy  [9]  in  conjunction  with  a  generalized  Mellin 
transform  to  analyze  stress  singularities  in  an  anisotropic  wedge.  It 
should  be  mentioned  that  plane  deformation  was  assumed  in  [7-9]  and  hence 
the  material  property  was  assumed  to  be  sy-mmetric  with  respect  to  the  plane 
of  deformation. 

In  this  paper  we  use  the  approach  of  Stroh  [10]  whose  analysis  was 
further  developed  by  Barnett  and  his  co-workers  (see  [11],  for  example)  to 
study  the  surface  waves  in  anisotropic  elastic  materials.  An  excellent 
review  article  on  surface  waves  in  anisotropic  elastic  materials  was  given 
by  Chadwick  and  Smith  [12] .  Although  no  stress  singularities  were  studied 
in  [10-12],  their  approach  is  used  here  to  find  the  stress  distribution  at 
an  anisotropic  wedge.  A  recent  study  by  Dempsey  and  Sinclair  [3]  on 
isotropic  elastic  wedge  problems  shows  that  the  singularity  analysis  can  be 
accomplished  without  resorting  to  the  Mellin  transform  even  when  the 
boundary  conditions  are  not  homogeneous  [4].  Following  their  analysis  and 
using  the  approach  of  Stroh,  we  present  here  possible  forms  of  stress 
distribution  near  the  vertex  of  a  wedge  or  a  composite  wedge  of  anisotropic 


materials. 
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The  stress  distribution  near  the  vertex  of  a  wedge  or  a  composite 
wedge  depends  on  whether  the  eigenvalues  p  of  the  elasticity  constants 
are  distinct.  It  also  depeods  on  whether  the  order  of  singularity  k  is 
a  single  or  multiple  root.  The  purpose  of  this  paper  is  to  show  how  one 
can  derive  the  form  of  stress  distribution  when  p  and/or  k  are  not 
single  roots.  We  also  show  the  geometrical  meaning  of  the  complex  variable 
2  in  terms  of  the  polar  coordinates  (r,9). 

Finally,  as  an  application,  we  consider  the  stress  singularity  at  a 
crack  tip  of  general  anisotropic  materials. 

2.  BASIC  EQUATIONS 

In  a  fixed  rectangular  coordinates  x^,  (i=  1,2,3),  let  u^, 

and  be  the  displacement,  stress  and  strain,  respectively.  The 

continuity  condition,  the  stress-strain  law  and  the  equations  of  equilibrium 
can  be  written  as 


c. .  =  (9u./9x.  +  du./dx.)/2 

ij  ^  i  j  j 

(1) 

'^ij  "  '^ijkil  ^kJl 

(2) 

9o. ./9x.  =  0 
ij  J 

(3) 

where 

'^ijkil  ”  ^jik2.  '  ^ijilk  ”  '^kilij 
are  the  elasticity  constants  of  the  anisotropic  material. 

otherwise,  repeated  indices  imply  summation. 

(4) 

Unless  stated 

We  assume  that  u^  and  are  independent  of  the 

Hence  we  assume  that 

Xj-coordinate. 

2  =  x^  ♦  px^ 

(5) 
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‘^ll 

"16 

"15 

Qik  = 

‘^ilkl 

"61 

"66 

"65 

l"51 

"56 

‘^55. 

"16 

"l2 

"14 

‘^ilk2 

"66 

"62 

"64 

"56 

"52 

"54. 

"66 

"62 

"64 

T  = 

“  ik 

‘^i2k2 

"26 

"22 

"24 

"46 

"42 

"44 

(10)  can 

then  be 

written  as 

D-u  = 
ik 

‘’ik  *  p"ik 

*  \i> 

P 

(14) 


(15) 


and  vanishing  of  the 

determinant  D.,  means 
ik 

"ll>2pci^*p2c66 

"16"P^"12*"66)"P^"26 

"l5"P^"l4""56^"P“"46 

"16"P^"12""66)"P^"26 

"66"2PC26-P^C22 

"56"P("25""46)"P^"24 

=0  (16) 

"l5*Pt"l4""56^"P^"46 

"56"P^"25*"46^*pS4 

"55"2P'^45*P^"44 

Eq.  (16)  provides  six  eigenvalues  of  p. 

For  each  of  p 

the  associated  u.'s 

1 

are  obtained  from  Eq 

.  (9).  In 

general,  (i  =  1 

2,3)  are  all  non- zero. 

Hence,  u^^ ,  u^  and 

u.  are 

0 

coupled. 

As  to  we  let  j  =1  and  2,  respectively,  in  Eq.  (8)  and  use 

the  notations  of  Eq.  (14)  to  obtain 

"ii  • 

^i2  '  ‘“ki  *  P'^k'^ 


C17) 
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It  follows  from  Eqs.  (9,15,17)  that 


and  hence 


"'■12  ^  ’^"^22  ■  ^11  P  ^22  ’  "'•13  ■^"'^23 


(18) 


(19a) 


Therefore,  of  the  six  components  jj  *  '*'®  J^eed  is  ’^22’  ^33  ^23‘ 

They  are  obtained  from  Eqs.  (17)  and  (8)  with  the  aid  of  Eq.  (14).  The 
results  can  be  casted  in  the  following  form: 

T.  =  (C.  ,  +  pc.  ,)u,  +  (c.  ,  +  pc.  .,)0-  +  (c.  .  +  pc.  ,)U,  , 

(i  =  2,3,4) 


where 


"'■2  ^  "'^22  ’  "^3  "‘'33  ’  "'■4  "'^23 


(19b) 

(19c) 


Notice  that  since  and  T^j^  are  symmetric,  so  is  D^j^.  .Notice 

also  that  c  ,  (j  =  1 ,2, . . . ,6)  are  not  present  in  Eq.  (16).  Therefore,  the 
eigenvalues  p  are  independent  of  these  elastic  constants.  In  fact,  the 
stress  singularities  are  also  independent  of  these  elastic  constants. 

Equation  (16)  is  a  sextic  equation  in  p.  If  the  strain  energy  is 
positive  definite,  it  can  be  shown  that  p  cannot  be  real  [10,12].  There¬ 
fore,  we  would  have  3  pairs  of  complex  conjugate  roots  for  p. 

Wien  the  material  property  is  symmetric  with  respect  to  the  (x, ,x,} 
plane  or  to  the  (Xj^.x^)  plane,  Eq.  (16)  reduces  to  a  cubic  in  p^  [19]. 
Since  every  cubic  has  at  least  one  real  root,  one  of  the  p's  will  be 
purely  imaginary  when  (x^.x^)  or  (Xj^,Xj)  is  a  plane  of  symmetry. 


4.  UNCOUPLING  OF  U3  FROM  AND  U2 

When  the  material  property  is  symmetric  with  respect  to  the  (x^^.x.,) 


plane,  we  have 


*^14  “  *^15  “  '^24  “  '^25  ^^34  “  '^35  “  *^46  "‘^56°° 


(20) 
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Equation  (16)  then  reduces  to 


Ci,,2pcj^.p'c^^ 

'l6"P'"l2*'=66>'P"‘26 

0 

'=66*"P<=26''’^‘=22 

0 

=  0 

(21) 

0 

0 

"55"2P<^45"P^^44 

Therefore,  instead  of  a  sextic  we  have  a  quartic  equation  and  a  quadratic 
equation  in  p. 


If  p  is  a  root  of  the  quartic  equation,  we  see  from  Eqs.  (9,21)  that 
0^=0.  Moreover,  Eqs.  (19)  show  that  “ ’'■23  "  ^ •  Hence,  we  have  a  plane 

deformation. 

Similarly,  if  p  is  a  root  of  the  quadratic  equation,  =  0  and 

^11  ”  ^22  ~  ^33  ~  ^12  ~  anti-plane  deformation. 

Therefore,  when  Eq.  (20)  holds,  the  plane  deformation  and  the  anti- 
plane  deformation  are  uncoupled.  Since  the  system  is  linear,  we  may  consider 
them  separately  when  Eq.  (20)  holds. 

5.  GEOMETRICAL  INTERPRETATION  OF  Z  =  + px2 

Let  a  and  B  be,  respectively,  the  real  and  imaginary  part  of  p  so 

that 

p  =  a  +  6i,  3>0  (22) 

We  assumed  B>0  because  the  conjugate  of  p  will  have  the  negative 
imaginary  part.  Using  the  polar  coordinates  with  the  origin  at  Xj^  =  X2  =  0, 
we  have 

x^  =  rcos9  ,  x^  =  r  sin9  .  (23) 

i>l^ 


Hence, 


Z  =  Xj^  +  px^  =  X  +  iY  =  rpe 


(24) 
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where  , 

X/r  =  COS0  +  cisinS  =  pcosijj 

Y/r  =  6sin9  =  psinij/ 

p^  =  (cose  +  asin0)^  +  g^sin^O 


(25) 


It  is  not  difficult  to  show  from  Eqs.  (23-25)  that  a  unit  circle  in  the 
(x^.x^)  plane  maps  an  ellipse  in  the  (X,Y)  plane.  Fig.  1.  If  the  (Xj^.x^) 
plane  is  a  stretchable  sheet,  one  obtains  the  ellipse  by  first  stretching 
the  circle  uniformly  3  units  in  the  x^-direction  and  then  shear  the 
sheet  with  the  x^-axis  fixed  until  point  b  displaces  a  unit  horizontally. 
From  point  a  in  (x^^.x^)  and  (X,Y)  planes  we  see  the  geometrical  relation¬ 
ship  between  9,  p  and  <p.  From  Eq.  (25),  notice  that  p  and  \p  depend  on 
9  and  p  but  are  independent  of  r.  Notice  also  that 

p=l,  ip  =  d  ,  at  9  =  0,  +TT  (26) 

If  p  is  purely  imaginary,  we  also  have,  In  addition  to  Eq.  (26), 

p  =  g,  =  0,  at  0  =  ±it/2,  ±3Tr/2  when  a  =  0  (27) 

For  isotropic  materials,  p=  +  i  is  a  multiple  root  of  order  3.  Thus 
the  ellipse  in  the  (X,Y)  plane  reduces  to  a  unit  circle.  Hence, 


and 


P  =  1  ,  =  9  , 

Z  =  Xj^  +  ix2  =  re^^ 


(28) 

(29) 


which  is  the  well-known  complex  coordinate  for  (Xj^,X2)  in  two-dimensional 
elasticity  problems  of  isotropic  materials. 


6.  STRESS  DISTRIBUTION  NEAR  THE  VERTEX  WHEN  p's  ARE  DISTINCT 

To  find  the  stress  distribution  and  the  stress  singularities  at  the 
vertex  of  a  wedge,  we  choose 
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f(2) 


_L 


(30) 


where  k  is  the  order  of  singularity  to  be  determined  by  the  boundary  condi- 
tions.  As  we  mentioned  earlier,  the  eigenvalues  p  are  all  complex  numbers 
and  consist  of  three  pairs  of  complex  conjugates.  In  this  section  we  assume 
that  the  eigenvalues  are  distinct.  Using  Eq.  (30)  in  Eqs.  (6,7)  for  all 
eigenvalues  and  forming  a  linear  combination  of  them  leads  to 


u. 

1 


a.  . 
ij 


(A^.fZ^' 


>  8^u.Z^‘‘")/(1-k:)  -  ... 

+  B.t.  +  . . . 

1  ij 


(31) 


^32) 


where  Aj^,B^,...  are  constants  which  may  be  complex  and  an  overbar  denotes 
a  complex  conjugate.  For  simplicity  only  the  terms  associated  with  one  pair 
of  eigenvalues  are  written  explicitly  to  avoid  introducing  an  additional 
subscript  for  the  eigenvalues.  The  dots  denote  terms  associated  with  the 
remaining  two  pairs  of  eigenvalues. 

It  should  be  pointed  out  that  Uj^  as  given  by  Eq.  (9)  is  not  unique 
and  can  have  an  arbitrary  multiplicative  constant.  The  constants  A^^  and 
B.  an  (31.32;  represent  this  arbitrary  multiplicative  constant. 

For  a  wedge  or  a  composite  wedge,  by  substituting  Eqs.  (51,52)  in  the 
boundary  conditions  (which  include  the  interface  conditions  if  the  wedge  is 
a  composite),  one  obtains  a  system  of  linear  algebraic  equations  in  A^  , 

B^  ,  .  .  . ,  which  may  be  written  as 


K. .c.  =  q. 
ij  j 


(33) 


where  is  a  square  matrix  which  depends  on  <,  c^  is  a  column  matrix 

whose  elements  are  Aj^.B^,...,  and  q^^  is  a  column  matrix  which  depends  on 
the  boundary  conditions.  If  the  boundary  conditions  are  homogeneous,  q.  =0. 
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In  this  case,  a  nontrivial  solution  exists  if  the  determinant  of  K. . 

ij 

vanishes.  The  roots  of  this  determinant  provides  the  values  for  k.  Let 


K  =  C  +  ni  ‘  (34) 

where  C  and  n  are  real.  If  0<5<1,  we  have  a  singularity  at  r  =  0. 

Since  u^  and  are  real,  only  the  real  parts  or  the  imaginary 

parts  on  the  right-hand  sides  of  Eqs.  (31,32)  should  be  considered.  They 
will  have  different  expressions  depending  on  if  the  root  <  is  real  or 
complex. 


A.  K  =  Real 


Since  and  are  in  general  complex,  let 


0  .  =  V .  e 
1  1 


xai 


T.  .  =  t.  .e 
ij  ij 


ib 


(35) 


where  v^,  a^,  t^^  and  are  real  and  repeated  indices  do  not  imply 
summation  here.  With  Eqs.  (24,35),  the  real  parts  of  Eqs.  (31,32)  can  be 
written  as 


=  (rp)^'^v^{M^cos[a^+(l-C)ip]+N^sin[a^+(l-5)i|^]}/(l-C)  +  .  .  .  (36) 

=  (rp)'^t^j{M^cos(b^j-5ip)+N^sin(b^^-C'^)}+.  •  •  (37) 


where  are  related  to  and  are  real.  The  imaginary 

parts  of  Eqs.  (31,32)  provide  no  new  expressions. 


B. 


assuming 


K  =  g  +  in,  complex 

When  K  is  a  complex  root  there  is  no  loss  in  generality  in 
n>0  because  if  k  is  a  root,  so  is  <.  We  then  have 


Z 


-K 


(38) 
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The  real  parts,  or  the  imaginary  parts  of  Eq.  (32)  now  become 

'^ij  ”  j  ^ e'^'^(Mj^cos(J)^j +N^sin(J)^j ) +e  ^'*^(M^cos<J)^j +N^sin(J)^j ) }+ .  . .  (39) 


where 


“^ij  b^j-C';'±nln(rp) 


(40) 


+ 

and  M",  are  real  constants  and  are  related  to  and  .  A  similar 

equation  may  be  written  for  u,  .  We  see  that  o.  .  is  oscilatorv  and 
unbounded  as  r-*-0.  As  expected,  Eq.  (39)  reduces  to  Eq.  (37)  when  p  =  0. 

In  the  sequel,  we  will  consider  only  the  cases  in  which  k  is  a 
complex.  The  solution  for  a  real  <  is  deduced  by  letting  n  =  0. 


7.  STRESS  DISTRIBUTION  .NEAR  THE  VERTEX  VJHEN  p  IS  A  DOUBLE  ROOT 

When  p  is  a  double  root  of  Eq.  (16),  we  have  only  two  pairs  of 
distinct  eigenvalues  instead  of  three.  It  is  not  difficult  to  see  that,  if 


u.  =  u.Z^‘'^/(l-<) 


-K 


a.  .  =  T. .Z 
ij  iJ 

are  the  solutions  corresponding  to  the  double  root  p,  so  are 


1  d  r  .,1-K-, 

^i  =  d^^^i"  ^ 


(41) 

(42) 


-L-ulz^-^ 

l-K  1 


.  0.  Z-x.  j 


a.  .  =  T.  .Z~*^  -  KT.  .Z''^~^Xt 
IJ  IJ  IJ  2 

where  a  prime  stands  for  differentiation  with  respect  to  p.  Since 


x^  =  (Z-Z)/(2p), 


we  have 


l-K  1  ..  =,-K 


u.  =  (t-^-  ^  ’J.)Z  -  ^  'J-3Z 

1  l-K  1  2p  1  2p  1 


'  K  ^  <  =,-K-l 

a..  =  (t..  -  -r-T..)Z 

IJ  IJ  :p  IJ-'  2p  IJ 


(43) 

(44) 

(45) 

(46) 

(47) 


-  •  f 
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is  obtained  by  differentiating  Eq.  (9)  with  respect  to  p: 


D-u'-’,' 
ik  k 


^lk\ 


0 


(48) 


The  existence  of  a  solution  for  and  from  Eqs.  (9,48)  will  not  be 

/ 

discussed  here.  Likewise,  is  obtained  by  differentiating  Eq.  (8)  with 

respect  to  p. 

Notice  that  and  obtained  from  Eqs.  (9,48)  are  not  unique. 

They  both  have  a  multiplicative  constant,  say  and  A2,  respectively. 

-  fC 

The  Z  term  in  Eq.  (47)  together  with  the  same  term  with  p  replaced  by 
p  is  essentially  similar  to  Eq.  (32)  and  hence  would  yield  an  expression 
similar  to  Eq.  (39).  The  last  term  in  Eq.  (47)  will  produce  a  second 
independent  solution.  This  is  obtained  from 


0.  .  =  A_  T.  .  Z2  +  B-  —  T.  .2Z 
ij  2  2p  IJ  2  2p  IJ 


(49) 


by  taking  the  real  or  imaginary  parts  of  the  right-hand  side.  Therefore, 
when  p  is  a  double  root,  we  have  the  following  second  independent  solution 
for  in  addition  to  Eq.  (39): 


0^^  =  (rp)  ^t^j  {e'^'^[M2C0s((J)^^  -  2\p)  +  N2sin(^)^j  -  2’1<)  ] 

+  e"'^'^[M2COs(4>^j  -  2^)  +  N^sinfc))^^  -  2iJ^)  ]  (50) 

where  <))*^  are  defined  in  Eq.  (40)  and  M^,  N*  are  related  to  A^,  p 
and  K.  Equation  (50)  applies  to  the  case  when  ic  is  complex.  For  a  real 
K,  we  simply  let  n=0. 


8.  STRESS  DISTRIBUTION  NEAR  THE  VERTEX  tVHEN  p  IS  A  TRIPLE  ROOT 

For  isotropic  materials,  p  =  ±i  is  a  triple  root.  However,  since  u^ 
is  uncoupled  from  u^  and  u^  for  isotropic  materials,  p  is  actually  a 
double  root  when  we  consider  u^^  and  U2  only.  Hence  the  previous  section 
on  a  double  root  p  applies  to  isotropic  materials. 


We  have  not  seen  an  example  other  than  isotropic  materials  for  which 
p  is  a  triple  root.  If  there  is  one,  and  if  is  not  uncoupled  from  u^^ 

and  u^,  we  see  that  a  thisd  independent  solution  is 


u. 

1 


^ 


(51) 

(52) 


Following  a  similar  procedure  in  deriving  Eq.  (50),  the  real  expressions  for 
the  third  independent  solution  when  p  is  a  triple  root  can  be  obtained 
from  Eq.  (50)  with  2ip  replaced  by  and  the  subscripts  2  replaced  by  3. 


9.  STRESS  NEAR  THE  VERTEX  WHE.V  K  IS  A  DOUBLE  ROOT 

Up  to  now,  we  tacitly  assumed  that  K  is  a  single  root  of  the 
determinant  of  and  hence,  other  than  a  multiplicative  constant,  the 

homogeneous  equation  of  Eq.  (33)  has  a  unique  solution  for  c^  whose 
elements  are  the  coefficients  If  <  is  a  multiple  root,  then 

may  not  be  unique  and  we  have  other  new  solutions. 

Let  K  be  a  double  root  of  the  determinant  defined  in  Eq.  (33) 

with  q^  =  0.  Then,  not  only  Eqs.  (31,32)  are  the  solutions,  but  also  are 


Since 


d  r  A  ,1-Ki  -  d  r  B  =1-<1 

u.  =  u.  2  }  +  u.  ^ 

1  X  dK  l-<  1  dtc  l-< 


a.  .  =  T.  .  ^{A  Z"^}  +  f .  .  ;^{B  +  . . . 

ij  ij  d<  ij  dK 


T.  .  ;^{A  Z*''}  =  ^  T.  .Z’*^  -  A  T.  .Z'^lnZ, 
ij  dK  dK  ij  ij 


the  first  term  on  the  right  is  essentially  the  same  as  the  first  term  of 


Eq.  (32).  The  second  term  provides  a  new  solution  for  a.  .  when  k  is  a 

ij 

double  root: 
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a  =  A-T.  .Z’'^ln  Z  +  B-T.  .2"'^ln  Z 

ij  2  ij  2  ij 


(56) 


The  real  or  imaginary  parts  of  Eq.  (56)  have  different  expressions  depending 
on  whether  p  is  a  single  root  or  a  multiple  root. 

A.  p  is  a  single  root 

When  p  is  a  single  root,  the  real  or  imaginary  parts  of  Eq.  (56) 
have  the  expression: 

=  (rp)“^t^j{e'^^[M^(ln(rD)cos(J)^T  -  ij;sin<()^T) 

+  N2(ln(rp)sin())^T  +  li^cost^)^?)  ] 

+  e’'^'^[M2(ln(rp)cos(J)^*  -  ij;sin())^p 

+  N2(ln(rp)sin({)^j  +  4)cos(})^p]}  (57) 

As  before,  4)*^  are  defined  in  Eq.  (40)  and  M;,  N*  are  related  to  A, 
and 


B.  p  is  a  multiple  root 

Let  us  consider  first  the  case  in  which  p  is  a  double  root. 
Then,  in  addition  to  Eq.  (56),  we  also  have  the  solution 


a..  =  A-  :i(T.  .Z''^lnZ)  +  :^(T.  •f'^lnZ) 

11  2  dp'-  ij  2  dp'-  ij 


However,  since 


^(T.  •Z^'^lnZ)  =  (T.'.  -  T.  .)Z‘'^ln2 
dp^  ij  ^  ^  ij  2p  ly 


+  T.  .  (Z  -  ZZ  )  +  t.  .ZZ  InZ 
2p  ir  2p  ij 


(58) 


(59) 


where  use  has  been  made  of  Eq.  (45),  only  the  last  term  provides  a  new  solu¬ 
tion.  The  rest  of  the  terms  in  Eq.  (59)  have  appeared  in  Eqs.  (56,32,47). 
Therefore,  a  new  solution  when  p  is  a  double  root  is 
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a.  .  =  A.,  ^  T,  .2Z''^"^ln  Z  +  B,  —  f.  .2Z''^'^ln  Z  (60] 

ij  2  2p  ij  2  2p  ij 

The  real  or  imaginary  parts^of  Eq.  (60)  have  the  expression  which  is  obtained 
from  Eq.  (57)  with  replaced  by  (tj)^^  -2’^). 

Similarly,  if  p  is  a  triple  root,  it  is  not  difficult  to  show  that 
the  new  solution  is  obtained  from  Eq.  (57)  by  replacing  by  (0:^  -  4i]/) . 

We  see  from  Eq.  (57)  that  has  the  singularity  of  r"^  In  r.  The 

existence  of  a  solution  of  Eq.  (57)  depends  on  the  existence  of  a  solution 
for  A  and  d.i/d<  in  Eq.  (55).  Since  A  is  an  element  of  c^  in 
Eq.  (55),  the  existence  of  A  and  dA/dx  depends  on  the  existence  of  a 
solution  for  c^  and  dc^/d<  from  the  following  equations 

K. .c.  =0  C61) 

IJ  J 

K. .(dc./d<)  +  (dK../dK)c.  =  0  (62) 

ij  J  ij  J 

A  discussion  of  the  solution  of  Eqs.  (61,  62)  can  be  found  in  [3]. 

10.  STRESS  NEAR  THE  VERTEX  WEN  K  IS  A  TRIPLE  ROOT 

When  <  is  a  triple  root,  one  can  follow  the  same  reasoning  as  in  the 
previous  section  for  a  double  root  k.  Therefore,  the  new  solution  for  a 
triple  root  <  is  obtained  by  replacing  d/dx  by  d^/dK^  in  Eqs.  (53,54). 
Equation  (56)  then  is  replaced  by 

°ij  =  A2T.jZ'^(lnZ)"  .Z'^^dn!)^  (63) 

and  Eq.  (57)  is  modified  by  replacing  In  (rp)  by  [In  (rp)  ]  ^  ,  which  is 

the  real  part  of  (InZ)^,  and  ip  by  2ij;  In  (rp),  which  is  the  imaginary  part 
of  (InZ)^ 
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11.  GENERAL  EXPRESSION 

We  can  sununarize  the  results  obtained  so  far  in  the  following  form. 

Let  n  be  the  multiplicity  of  p  and  m  be  the  multiplicity  of  k.  If 
P  ^ 

we  write 

a..  =  r'^F^.(r.8) 


then  F  consists  of  a  linear  combination  of  the  real  and  imaginary  parts 

ij 

of  the  following  expression 


t..0'V’^^{lnCrp)±i^l'"'^{cos[b..-^4^-2Cn-l)'l-^nlnCrp)] 

+  isin  [b^.-Cii)-2(n-i)ii/^nln(rp)]} 

for  each  p  and  for  all  integers  m  and  n  subjected  to  the  limitations 

As  we  stated  before,  P  and  ijj  depend  on  0  but  not  on  r. 


[64b) 


(64c) 


12.  SINGULARITY  AT  A  CRACK  TIP  FOR  ANISOTROPIC  SOLIDS 

Consider  an  infinite  anisotropic  solid  with  a  crack  plane  which  is 
located  at  <  0  of  the  (x^.xp  plane.  Hence,  ^2^=0.  [j  =1,2,3)  at 
0  =  ±7r.  Using  Eq.  (26),  Eq.  (37)  for  9  =  tt  and  -tt  reduces  to 


t2j{MiCOs(b.,j-CTr)  +  NjSin(b2j -Ci^) }  +  •••  =0 

t2.{M^cos(b2.+CiT)  +  N^sin(b2j+5^)}  +  •••  =0 
^  ^  (j  =  1.2,3) 


If  we  set  5  =  %,  we  have 

t2jlM^sin(b2j)-N^cos(b2j)}  +  ••*  =0 

t2.{-MiSin(b2j)+NiCOs(b2j)}  +  •••  =0 

(j  =  1.2.3) 


(65) 


(66) 
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Equation  (66)  consists  of  6  equations  for  and  can  be  written  in 

the  form  of  Eq.  (33)  with  =  0.  Since  the  first  three  equations  are 
identical  to  the  last  three  equations,  C  =  *5  is  a  triple  root  of  the  deter¬ 
minant.  We  can  therefore  let  5  =  4.  H  =  0,  =  3  in  Eqs.  (64).  Disregarding 

the  dependence  on  0,  the  singularities  at  the  crack  tip  in  a  general 

anisotropic  material  are  r'^  and  possibly  r~^  In  r  and  r“'^(lnr)^.  The 
-i'  -i' 

existence  of  r  *lnr  and  r'*(lnr)^  depends  on  the  existence  of  a 
solution  for  c^ ,  dc^/dK,  d^c^/d<^  from  Eqs.  (61,62)  and  an  equation 
obtained  by  differentiating  Eq.  (62)  with  k. 

13.  SUMMARY  A-MD  CONCLUSION 

We  have  presented  here  a  means  to  determine  the  order  of  singularity  < 
at  the  free-edge  of  an  interface  in  a  layered  composite  in  which  each  layer 
is  anisotropic.  Although  the  order  of  singularity  does  not  depend  on  the 
stacking  sequence  of  the  layers  in  the  composite,  the  coefficients  of  the 
singular  terms  which  are  related  to  the  intensity  factor  do.  These  co- 
efficients  can  be  determined  only  by  solving  the  complete  boundary  value 
problem.  One  may  use  a  special  finite  element  at  the  free-edge  using  the 
analyses  presented  here  and  regular  finite  elements  elsewhere  in  solving  the 
complete  boundary- value  problem. 
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CHAPTER  6 

TRANSIENT  WAVE  PROPAGATION  NORMAL 
TO  THE  LAYERING  OF  A  FINITE  LAYERED  MEDIUM 


ABSTRACT 

Plane  wave  propagation  in  the  direction  normal  to  the  layering  of  a 
periodically  layered  medium  is  studied.  A  period  consists  of  two  layers 
of  homogeneous,  isotropic,  linear  elastic  or  viscoelastic  materials.  The 
layered  medium  is  of  finite  extent  and  hence  consists  of  a  finite  number 
of  layers.  A  theory  is  presented  by  which  the  layered  medium  is  replaced 
by  an  "equivalent"  linear  homogeneous  viscoelastic  material  such  that  the 
stress  or  the  velocity  in  the  latter  and  in  the  layered  medium  are  identical 
at  the  centers  of  the  alternate  layers.  Transient  waves  in  the  layered 
medium  are  then  obtained  by  solving  the  transient  waves  in  the  "equivalent" 
homogeneous  viscoelastic  medium.  Solutions  at  points  other  than  the  centers 
of  the  alternate  layers  are  also  presented.  Numerical  examples  are  given 
for  transient  waves  in  an  elastic  layered  medium  due  to  a  step  load  applied 
at  one  of  the  boundary  while  the  other  boundary  is  fixed.  Comparisons  with 
the  exact  solutions  by  the  ray  theory  show  that  the  present  theory  can  predict 
very  satisfactorily  transient  waves  in  a  finite  layered  medium.  The  theory 
of  viscoelastic  analogy  applies  to  other  cases  for  which  exact  solutions  by  the 
ray  theory  are  not  available,  such  as  the  case  of  finite  layered  medium 
with  prescribed  boundary  conditions  which  are  time -dependent . 
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1.  INTRODUCTION 

Most  of  the  approximate  theories  for  wave  propagation  in  a  layered 
medium  focus  on  the  determination  of  the  dispersion  relation  or  the  fre¬ 
quency  equation  due  to  a  harmonic  oscillation  [1-4],  although  some  of  the 
theories  are  able  to  predict  the  late-time  asymptotic  solution  in  a  semi¬ 
infinite  layered  medium  due  to  a  step  load  applied  at  the  boundary.  For 
the  latter,  exact  theories  may  be  used  to  find  the  asymptotic  solution  and 
the  wave-front  solution  [5-7]  . 

To  predict  the  transient  response  at  points  not  necessarily  far  away 
from  the  impact  end  (where  the  asymptotic  solution  does  not  apply)  and  to 
points  not  necessarily  near  the  wave-front,  a  new  theory  based  on  the 
analogy  between  the  dynamic  response  of  a  semi-infinite  layered  medium  and 
a  semi-infinite  homogeneous  viscoelastic  medium  has  been  proposed  recently 
by  Ting  and  Mukunoki  [8].  The  fundamental  idea  is  to  characterize  the 
layered  medium  by  an  "equivalent"  homogeneous  viscoelastic  medium  such  that 
the  dynamic  response  of  the  latter  is  identical  to  that  of  the  layered  medium 
at  the  centers  of  the  alternate  layers.  Although  the  idea  of  modeling  a 
composite  by  a  viscoelastic  medium  is  not  new  [9,10],  the  "theory  of  visco¬ 
elastic  analogy"  introduced  in  [8]  succeeds  in  obtaining  the  exact  form  of 
the  "equivalent"  relaxation  function  for  the  layered  medium. 

Since  wave  propagation  in  a  homogeneous  linear  viscoelastic  medium  can  be 
solved  easily  by  many  known  numerical  schemes  (see  [11],  for  e.xample),  one 
can  obtain  the  transient  wave  solution  in  a  layered  medium  by  solving  the 
transient  waves  in  the  "equivalent"  homogeneous  viscoelastic  medium. 

The  layered  medium  considered  in  [8]  is  of  semi-infinite  extent.  In 
this  paper  we  extend  the  theory  to  the  case  of  a  finite  layered  medium. 

First,  we  derive  the  general  solution  in  the  form  of  Laplace  transform  for 
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waves  propagating  normal  to  the  layerings  of  a  finite  layered  medium.  The 
general  solution,  which  is  applicable  to  any  point  in  the  finite  layered 
medium,  contains  two  arbitrary  coefficients  which  can  be  determined  from 
the  boundary  conditions  of  the  finite  layered  medium.  Next,  we  apply  the 
general  solution  to  certain  points  in  the  layered  medium,  namely,  the  centers 
of  each  layer.  We  show  that  the  general  solution  at  the  centers  of  each 
layer  is  analogous  to  the  general  solution  for  waves  propagating  in  a  homo¬ 
geneous  viscoelastic  medium.  From  this  analogy  we  obtain  the  viscoelastic 
relaxation  function  of  the  "equivalent"  homogeneous  viscoelastic  medium. 
Several  analogies  can  be  made  depending  on  whether  one  is  interested  in  the 
stress  response  or  the  velocity  response  in  the  layered  medium.  The 
analogies  obtained  here  are  more  general  than  that  presented  in  [8]  and 
can  be  applied  to  the  semi-infinite  medium  as  well.  In  finding  a  means 
for  determining  the  response  at  points  other  than  the  centers  of  the  layers, 
we  inadvertently  obtain  a  characteristic  relation  in  an  integral  form  for 
one-dimensional  waves  in  homogeneous  viscoelastic  media.  In  the  literature, 
this  is  in  a  differential  form. 

2.  BASIC  EQUATIONS 

Consider  a  periodic  layered  medium  as  shown  in  Fig.  1  in  which  each 
period  2a)  consists  of  two  layers  of  homogeneous,  isotropic,  linear  elastic 
or  viscoelastic  materials.  The  two  different  materials  in  the  layers  will 
be  designated  as  material  1  and  2,  respectively.  Thus  material  1  occupies 
layers  1,  3,  5...  while  material  2  occupies  layers  2,  4,  6...  The  thick¬ 
nesses  of  individual  layers  are  denoted  by  2h^  [i=l,2)  where  the  subscripts 
1  and  2  refer  to  material  1  and  2,  respectively.  We  will  assume  that  the 
layered  medium  is  initially  at  rest  and  occupies  the  region  0  <  x  <  1.  We 
choose  the  central  surface  of  layer  1  as  x  =  0  and  the  other  boundary,  x=  . , 
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is  assumed  to  be  at  the  central  surface  of  layer  N  where  N’  can  be  an  even 
or  odd  integer.  Hence, 

^  =  (N  -  1)0)  CD 

IVe  will  consider  plane  wave  propagation  in  the  direction  x  in  which  the 
only  non-vanishing  component  of  the  displacement  is  in  the  x  direction. 

We  therefore  have  a  one-dimensional  wave  propagation  problem  in  which  the 
equation  of  motion  and  the  continuity  of  the  displacement  are  given  by 

9a 

.  (i  =  1.2)  (2) 

9v. 

(i-1.2)  (5) 

where  a  dot  stands  for  differentiation  with  respect  to  the  time  t,  and 
"^i’  ^i’  '^i’  *^i  normal  stress,  normal  strain,  particle 

velocity  and  mass  density,  respectively.  Let  X^(t)  and  y^^Ct)  be  the  re- 
la,\ation  functions  of  the  materials.  For  elastic  materials,  X^(t)  and 
y^(t)  are  independent  of  t  and  are  identified  as  Lame  constants.  The 
stress-strain  relation  can  be  written  in  the  form  of  Stieltjes  convolution 

P^Cx.t)  =  J  g^(t  -  t')  d£.  (x,t')  ,  (4) 

0' 

g^Ct)  =  X.Ct)  >  2u.Ct)  ,  (5) 

where  we  have  assumed  that 

a.(x,0‘)  =  v.(x,0")  =  e^(x,0')  =  0  .  (6) 

3.  GENERAL  SOLUTIONS 

The  general  solution  to  Eqs.  (2-6)  can  be  obtained  by  the  method  of 
Laplace  transform  and  by  the  use  of  the  Floquet  theory.  We  define  the  Laplace 
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transfonn,  f(p),  of  a  function  f(t)  by 
f(p) 


■i: 


fCt)e‘P‘^  dt 


C7) 


After  applying  the  Laplace  transform  to  Eqs.  (2-6),  the  general  solution  for 
the  stress  and  the  velocity  in  layers  1  and  2  can  be  written  as 


5^(x,p) 


V^Cx.p) 


V2(x,p) 


cosh  (kjX)  +  sinh  (k^x) 

—  /  a,  sinh  (k.x)  +  B,  cosh  (k.x)}- 
^  1  1  1  1  ^ 

A^  cosh  (k^x  -  k^oi)  +  sinh  (k^x  -  k^co) 

^  |a2  sinh  (k^x  -  k^co)  +  B^  cosh  [k^x  -  k^u)! 


(8a) 

(8b) 

(8c) 

(8d) 


where 


0)  =  h^  +  h2 


k.  = 


(9) 


""i  “  '^i  •'Pi  PSi 


A^  and  (i  =  1,2)  are  determined  by  the  continuity  condition  at  x  =  hj^ 


r  - 

(h^,p)  = 

^2 

■  • 

> 

_ 1 

(h^,p) 


(10) 


and  the  quasi-periodicity  property  of  the  solution  together  with  the  con¬ 


tinuity  condition  at  x  =  2a)  -  h 


1 


'^2‘ 

(2a)-h^,p)  = 

S' 2  . 

."1. 

(-hj,p) 


-2a)K 


(11) 


where  k  is  the  characteristic  exponent  [12].  Substitution  of  Eqs,  (8)  into 
Eqs.  (10)  and  (11)  leads  to  four  homogeneous  equations  for  A^  and  B^.  The 
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requirement  for  a  non-trivial  solution  results  in  the  following  equation 
for  the  characteristic  exponent  <: 


cosh(2a)K)  =  9  cosh(2kj^hj^  +  2k2h2)  -  (9-1)  cosh(2kj^hj^  -  2k2h2) 


Q  1/1  T  1 

9=t-(  —  +  2+  — 
4  \ 


Moreover,  A.  and  B,  are  related  bv 
11 


^  =  pMe'^ 

•s 


B 


where 


1 


—  =  -pLi  , 


=  m^^RpM 


^2  , 

—  =  -pL,  . 

^2 


L2  =  m2R/(pM) 


pM  = 


'"l^l^2  *  "'2^1^2 
cosh  (ojK) 


m2  cosh  (ojk) 
m2^j^^2  ^  m2^S^S2 


pR  = 


'"l*^1^2  "  '"2'^2^1 
m^m2  sinh(aK) 


sinhCcjK) 


1^2  "■’"1^2^  1 


=  cosh(k^h^)  ,  =  sinh(k^h^) 


(12) 

(13) 


(14) 


(15) 


Notice  that  if  we  interchange  the  subscripts  1  and  2,  the  expression  for  pR 
remains  unchanged  while  pM  becomes  (pM)  Therefore,  we  can  obtain  the 
Stieltjes  inversion  of  M(t)  by  simply  interchanging  the  subscripts  1  and  2 
in  the  expression  for  pM  and  applying  the  Laplace  inverse  transform. 

With  Eq.  (14) ,  the  general  solution  in  the  layers  1  and  2  as  expressed 
by  Eq.  (8)  can  now  be  reduced  to  a  solution  containing  only  one  coefficient, 
say  Aj^.  The  solutions  in  other  layers  are  obtained  by  the  quasi-periodicity 


relation: 


18S 


-2na»: 


C2na)+  x.p)  =  _  (x,p)  e 


where  n  is  an  integer.  Moreover,  we  see  from  Eq.  (12)  that  if  k  is  a 
characteristic  exponent,  so  is  -k.  Therefore,  in  addition  to  the  general 
solution  with  as  the  coefficient,  we  obtain  the  second  general  solu¬ 
tion  by  changing  the  sign  of  K.  The  coefficient  of  this  second  solution 
will  be  denoted  by  .  Consequently,  the  general  solution  for  the  stress 
and  velocity  at  any  point  x  in  the  layered  medium  can  be  written  as,  using 


Eqs.  (8,14,16), 


“  “  (  “  ^  ’*2n0i)K 

(2na)  +  Xj^ ,p)  =  < cosh(kj^Xj^)  -  pL^^  sinh(kj^Xj^)  >  e 

•f  Aj  I  cosh(kjX^)  +  pLj  sinh(k^x^)  (17a) 


where 


v^(2n(u  +  x^,p)  =  jjp-|sinh(kj^Xj)  -  pL^  cosh(k^x^)  |  ( 
A* 

+  J-|sinh(k^x^)  +  pL^  cosh(k^x^)|( 


-2ncoic 


(17b) 


a2(2na)  +  00  +  X2,p)  =  Aj^pM  |cosh(k2X2)  -  PL2  sinh(k2X2)  | 

+  A^pM  |cosh(k2X2)  +PL2  sinh(k2X,) (17c) 


V2(2na)  +  00  +  X2,p)  =  jjp  pM  |sinh(k2X2)  -  PL2  cosh(k2X2)  le' 


(2n+l)ajK 


A 

~  pM  |sinh(k2X2)  +  PL2  cosh(k2X2)  (17d) 


-  h^  <  x^  <  h^  ,  (i  =  1,2) 
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lifhen  proper  values  for  n  and  (or  x^)  are  chosen,  Eqs.  (17)  can  be  used 
to  determine  solution  at  any  point  in  the  layered  medium.  The  two  co¬ 
efficients  and  are  determined  from  the  boundary  conditions  at  x  -  0 
and  X  =  Z. 

In  the  next  section  we  will  show  how  one  can  obtain  the  solution  at 
the  centers  of  the  layers  by  solving  the  wave  propagation  problem  in  an 
"equivalent"  homogeneous  viscoelastic  medium.  Having  found  the  viscoelastic 
analogy  for  the  solution  at  the  centers  of  the  layers,  we  then  show  how  one 
can  obtain  the  solution  at  points  other  than  the  centers  of  the  layers  in 
terms  of  the  solution  at  the  centers  of  the  layers. 


4.  SOLUTION  AT  CENTERS  OF  LAYERS 


The  stress  and  velocity  at  the  centers  of  the  layers  have  specially 
simple  forms.  By  letting  x^  =  X2  =  0  in  Eq.  (17),  we  have 

Oj  (2n(i),  p)  =  A^e  +  A^e 

-  ,  ^^1  (  X  -2nax  r»  2nax\ 

v^(2na),  p)  =  —  (-A^e  .  A^e  ) 

52(2nu).<.,p)  =  pM(A^e-^-""l^<^  . 

pM(-Aje-(-""^^‘^  *  A;e^2n.l)uK^ 


(19) 


P^2 

v,(2nu)  +  0),  p)  = 
l  ra2 


We  now  consider  a  homogeneous,  isotropic,  linear  viscoelastic  medium 
which  occupies  0 < x < £  and  which  is  at  rest  at  t  =  O"  and  is  subjected 
to  certain  prescribed  boundary  conditions  at  x  =  0  and  x  =  Jl.  Let  $,  n 
and  V  be  the  normal  stress,  normal  strain  and  particle  velocity,  respec¬ 
tively.  Also,  let  p  and  G  be  the  "equivalent"  mass  density  and  the 
"equivalent"  relaxation  function  of  this  homogeneous  viscoelastic  material. 


•  t 

t 
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The  equation  of  motion,  the  continuity  condition 
and  the  initial  conditions  are 


a* 


3x 


n 


$(x,t) 


<&(x,0') 


l.t 

G(t  -  t')  dn  (x.t') 

0" 

=  v(x,o’)  =  n(x,o")  =  0 


,  the  stress-strain  relation 


(20) 


By  applying  the  Laplace  transform  to  Eqs.  (20),  the  general  solution  for  the 
stress  and  velocity  will  contain  the  exponential  term 

exp  /pp/G  (21) 

In  view  of  the  exponential  terms  in  Eqs.  (19),  we  will  define  the  "equivalent" 
relaxation  function  G(t)  by  the  relation 

<  =  /pp/G  (22) 

We  will  also  define  the  "equivalent"  mass  density  p  by  the  average  mass 
density  in  the  layered  medium  [4,8]; 

p  =  (p^h^  +  P2h2)/(h^+  h^)  (25) 

With  Eq.  (22),  the  general  solution  to  Eq.  (20)  can  be  written  as 

^(x,p)  =  a  e  +  a  e  (24a) 

V(x,p)  =  ^  a  +  i’  e'^'')  (24b) 

where  a  and  a'  are  arbitrary  functions  of  p. 

There  are  several  ways  to  identify  the  analogy  between  Eqs.  (19)  and  (24). 
If  the  stress  in  material  1  is  of  main  interest,  we  may  set 


lyi 


Ai  =  a 


r  I  -  I 

Aj^  =  a 


(25) 


we  then  have 


and 


a^(x,p)  =  $(x,p) 
v^(x,p)  =  pj^V(x,p) 

a^(x,p)  =  pMl(x,p) 
v^Cx.p)  =  pM  {pJ^VCx.p)} 


for  X  =  2na» 


(26a) 


for  X  =  (2n+l)uj 


(26b) 


where 


ppL. 

J.  =  - - 

1  Km. 


(i  =  1.2) 


(27) 


It  should  be  pointed  out  that  while  $  and  V  as  given  by  Eq.  (24)  are 
defined  for  all  x,  Eqs.  (26a)  and  (26b)  apply  only  to  x  =  2nu  and 
X  =  (2n+l)a),  respectively.  By  using  the  identity, 

j^/j^  =  (pM)^  =  012^/(111^12) 


the  last  of  Eq.  (26b)  can  be  written  as 

v^(x,p)  =  {p>j,V(x,p)}  ,  X  =  (2n+l) 

oM 


(28) 


u 


(29) 


With  Eq.  (29),  we  rewrite  Eqs.  (26)  in  the  following  form: 


aj^(2nu,t)  =  <f(2n(i),t) 


v^(2n(i),t)  =  V*(2na),t) 


02  (2nu)  +  u»,t)  = 


M(t  -  t')  d<t  (2n«  +  w,t') 


V2(2n(j  +  a),t)  =  f  M"^  (t  -  t ' )  dV*  (2nu)  +  u, t ' ) 


(50a) 

(30b) 

(30c) 

(30d) 


where 


V  (x,t)  = 


^0- 


J^(t  -  t’)  dV  (x,t') 


(30e) 
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and  M”^  is  the  Stieltjes  inverse  of  M.  (See  the  discussion  following  Eq.  (15) 
regarding  the  Stieltjes  inverse  of  M.)  Thus  the  stress  and  velocity  at  the 
centers  of  the  layers  are  related  to  the  stress  $  and  velocity  V  in  the 
"equivalent”  homogeneous  medium.  In  particular,  the  stress  at  the  centers 
of  the  odd  layers,  a^(2u,t),  is  identical  to  the  stress  in  the "equivalent" 
homogeneous  viscoelastic  medium. 

As  an  illustration  for  the  theory  of  viscoelastic  analogy,  we  consider 
an  elastic  layered  medium  which  is  fixed  at  the  center  of  the  14th  layer 
(i.e.,  2,  =  13u))  and  subjected  to  a  unit  step  stress  applied  at  x  =  0,  Since 
the  14th  layer  is  occupied  by  material  2,  we  have  the  following  boundary 
conditions ; 

0^(0, t)  =  H(t)  ,  v^ClSoj.t)  =  0  (31a) 

where  H(t)  is  the  Heaviside  step  function.  In  view  of  Eqs.  (50),  the  corre¬ 
sponding  boundary  conditions  for  the  "equivalent"  viscoelastic  medium  are: 

<f(0,t)  =  H(t)  ,  V  (15oj,t)  =  0  (31b) 

We  now  replace  the  elastic  layered  medium  by  the  "equivalent"  homogeneous 
viscoelastic  medium  whose  mass  density  p  and  the  relaxation  function  G(t) 
are  given  by  Eqs,  (23)  and  (22).  Because  of  the  complicated  expression  for 
K  as  given  by  Eq.  (12),  analytical  inversion  of  the  Laplace  transform 
G(p)  from  Eq.  (22)  does  not  appear  feasible.  We  therefore  resort  to  a 
numerical  Laplace  inversion  of  G(p),  [13].  The  result  is  shown  in  Fig.  2 
along  with  the  physical  parameters  of  the  elastic  layered  medium  used  in  the 
calculation.  The  physical  parameters  are  taken  from  [4].  Unlike  for  most 
real  viscoelastic  materials,  the  relaxation  function  for  the  "equivalent" 
viscoelastic  medium  is  not  a  monotonically  decreasing  function  of  t.  This 
was  also  predicted  by  Christensen  [10]  based  on  the  dielectric  theory.  In 
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[8]  one  can  find  a  discussion  on  the  behavior  of  G(t)  as  t-^-O  and  as 

well  as  the  value  of  G(t)  at  t  =  0. 

With  G(t)  given  by  Fig.  2  and  the  boundary  conditions  given  by 
Eq.  C31b) ,  we  integrate  Eq.  (20)  numerically  by  the  method  of  characteris¬ 
tics  [14]  for  the  stress  $  and  velocity  V  in  the  "equivalent"  homogeneous 
viscoelastic  medium.  The  stress  and  velocity  at  the  centers  of  the  layers 
in  the  layered  medium  are  then  determined  by  using  the  viscoelastic  analogy 
Eqs.  (30).  In  Fig.  3  we  present  ^>(4a),t)  which  is  the  stress  history  at 
the  center  of  the  5th  layer.  For  this  example,  the  exact  solution  in  the 
elastic  layered  medium  using  the  ray  theory  can  be  obtained  numerically  by 
keeping  track  of  every  reflected  and  transmitted  waves  at  the  interfaces 
of  the  layers  [8].  This  exact  solution  is  also  shown  in  Fig.  3  for  compari¬ 
son.  It  is  seen  that  the  agreement  is  excellent. 

In  Fig.  3  we  also  show  the  solution  obtained  by  the  effective  modulus 
theory  [15].  With  this  theory,  the  elastic  layered  medium  is  replaced  by 
a  homogeneous  elastic  medium  whose  effective  modulus  is  given  by 


and  whose  effective  mass  density  is  identical  to  the  "equivalent"  mass  den¬ 
sity  defined  in  Eq.  (23) . 

In  Figs.  4  and  5  we  show,  respectively,  the  velocity  history  at  the 
center  of  the  5th  layer  and  the  stress  history  at  the  center  of  the  8th 
layer  by  using  Eqs.  (30b)  and  (30c).  Since  $(x,t)  and  V(x,t)  have  already 
been  determined,  all  we  need  is  the  functions  Jj^(t)  and  M(t)  which  are 
defined  in  Eqs.  (27)  and  (15).  Jj^(t)  and  M(t)  are  obtained  numerically  by 
inverting  their  Laplace  transforms.  Again,  the  solutions  by  the  ray  theory 


'j? 
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and  the  effective  modulus  theory  are  also  shown  in  the  figures  for  compari¬ 
son. 

The  function  as  well  as  functions  M,  R,  and  L^,  (i  =  1,2)  defined 
in  Eq.  (15)  are  called  the  "auxiliary"  functions.  Like  G(t),  the  auxiliary 
functions  depend  only  on  the  physical  properties  and  the  geometrical  layer¬ 
ing  of  the  layered  medium.  They  are  independent  of  the  boundary  conditions. 
For  the  viscoelastic  analogy  given  by  Eqs.  (30),  only  the  functions  M, 
and  M  ^  are  needed.  Of  course,  if  a^(2uj,t)  is  the  only  quantity  desired, 
no  auxiliary  functions  are  needed. 

Before  we  study  the  solution  at  points  other  than  the  centers  of  the 
layers,  we  will  discuss  other  forms  of  viscoelastic  analogy  in  the  next 
section. 

5.  OTHER  FORMS  OF  VISCOELASTIC  ANALOGY 

In  Eqs.  (30)  we  present  one  form  of  viscoelastic  analogy  between 
Eqs.  (19)  and  (24).  The  analogy,  Eqs.  (30),  is  convenient  for  the  case 
when  the  stress  in  material  1  is  of  main  interest  because  according  to 
Eq.  (30a)  is  identical  to  $.  If  the  stress  in  material  2  is  of  main 
interest,  then  the  analogy  given  by  Eq.  (30c)  requires  a  convolution  inte¬ 
gral  with  the  auxiliary  function  M(t) . 

There  are  of  course  other  forms  of  viscoelastic  analogy  which  would 
be  more  convenient  for  other  situations.  If  the  stress  in  material  2  is 
of  main  interest,  one  may  set 

pMA^  »  a 

PMa|  =  a' 

Then  the  analogy  between  Eqs.  (19)  and  (24)  can  be  written  as 
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a,  (2nw  +  oj,  t)  = 
V,  (2na)  +  oj.t)  = 


a^(2n(jj,t) 


•'0- 


$(2n(i)  +  cjJ.t) 

V*  (2n(jj  +  a),t) 

M  ^  (t  -  t')  d$  (2nu),t') 


(2na),t) 

V*(x,t)  = 


•t 

=  M(t  -  t')  dV*  (2nai,t') 
■^0“ 

.t 

J,(t  -  t'}  dV  (x,t') 

O' 


(34) 


With  this  analogy,  the  stress  i  in  the  "equivalent"  viscoelastic  medium  is 
identical  to  at  x  =  (2n+l)a). 

Likewise,  if  the  velocity  in  material  1  is  of  main  interest,  the  visco¬ 
elastic  analogy  can  be  written  as 


Vj^C2noj,t)  =  VC2naj,t} 
aj^(2nu),t)  =  $*(2n(jj,t) 


v^  (2na)  +  a),t) 


,t 

M  ^  (t  -  t ' )  dV  (2na)  +  co  ,  t ' ) 
■'O" 


a,  (2na)  +  ai,t) 
t 

$*(x,t}  = 

■^0' 


M(t  -  t' )  d'l'* (2nw  +  (j0,t' ) 


J‘^t  -t')  d$(x,t') 


/ 


(55) 


Finally,  if  the  velocity  in  material  2  is  of  main  interest,  we  can  write 


V2(2na)  +  a),t)  =  V(2na)  +  u,t) 

02  (2n(jj  +  aj,t)  =  $*  (2nu)  +  uj,t) 
t 

v^(2na),t)  =  I  M(t  -  t')  dV  (2na),t') 
0" 


(36) 
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y^(2n<jj,t)  =  M‘^(t  -  t')  d4>*  (2nu),t') 


> 


i  (x,t)  = 


J‘^(t  -  t')  d4>(x,t') 


(36) 
Cont '  d 


Sometimes  the  boundary  conditions  may  influence  the  choice  of  a  visco¬ 
elastic  analogy.  For  example,  suppose  that  one  is  interested  in  the  stress 
in  material  1  for  the  problem  in  which  the  velocity  is  prescribed  at  x  =  0 
(i.e.,  v^(0,t)  is  known),  and  the  other  boundary  x  =  i  is  fixed.  We  could 
use  either  the  analogy  Eqs.  (30)  or  the  analogy  Eqs.  (35).  If  we  use 
Eqs.  (30),  we  obtain  directly  from  ^  but  then  we  have  to  transform 
the  boundary  condition  v^(0,t)  to  V(0,t)  by  using  Eqs.  (30b, e), 


V(0,t) 


J'^(t  -  t')  dv^(O.t') 


(37) 


0" 


before  we  solve  for  <J>  and  V  in  the  "equivalent"  homogeneous  viscoelastic 
medium.  If  we  use  Eqs.  (55),  we  can  solve  for  'I'  and  V  immediately  since 
V(0,t)  =  v^(0,t),  but  to  obtain  0^  from  $  a  convolution  integral  is 
required. 


6.  SOLUTION  AT  ARBITRARY  POINTS 


If  we  solve  for  A^  and  Aj^  from  the  first  two  equations  of  Eqs.  (19) 
and  substitute  the  results  into  Eqs.  (17a, b),  we  have 

kiXi  -ki^i 


1  /  ^1*1  "*^l^l\  - 

a^(2n(jj  +  x^,p)  =  +  e  ja^(2na),p) 

1  {  ’'1^1  ■'^l^l\  - 

*  I\®  ■  ®  /  '"iV^(2nw,p) 

1  t  ’^1^1  -kiX,\  _ 

v^(2nu+Xj,p)  =  j\e  +  e  ^  ^yv^(2nu),p) 


k. x.  -k, X, 
+  -^  le  ^  -  e 


(38a) 


(58b) 
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Similar  results  can  be  obtained  for  a,  (2nu  +  oj  +  x-,  ,p)  and  (2nw  +  uj  +  x^  ,pj . 

We  will  define  the  Laplace  transform  of  the  functions  , 

Ej^(x^,t)  and  Fj^(x^,t]  by 


D^(x^,P3 


1  ■’"r’^i 
—  e  ^  ^ 


E^Cx^.p)  = 


Fj_Cx^,p) 


-k^Xj 

pm^ 


(59) 


Ej^  and  Fj^  have  Che  following  physical  interpretation.  Suppose  that 
material  1  occupies  the  semi-infinite  space  x  ^  0  and  is  initially  at  rest. 
Then  Dj^(x^,t)  and  Ej^(Xj^,t)  are,  respectively,  the  stress  and  velocity 
history  at  x  =  x^^  due  to  a  unit  step  normal  stress  applied  at  x  =  0. 
Fj^(Xj^,t)  is  the  stress  history  at  x  =  x^  due  to  a  unit  step  velocity 
applied  at  x  =  0.  We  now  rewrite  Eq.  (38a)  as 


where 


1  = 


(2nco  +  x^  ,p)  =  y  (2na>,p)  e 


'1 


px^/c 


10 


1 


P  D^(-^l,p)  e 


•px^/c 


10 


+  -  ja^  (2na),p)  e 


■px^/c 


10 


p jD^(x^,p)  e 


px^/c 


10 


1  - 


+  y ]v^(2nu),p)  e 


pxj/c 


10 


P  )  F^(-x^,p)  e 


-px^/c 


10 


1 


-  V  (2nco,p)  e  ‘ 


^-pxf/c 


P  jF^Cx^.p)  e 


px^/Cio 


Cio  =  ^  gi(0)/p.  ,  (i  =  1.2) 


C40) 


(41) 


and  Dj^(-Xj^,t)  is  obtained  from  Dj^(x^,t)  by  analytically  extrapolating  the 
later  from  Xj^  >  0  to  Xj^  <  0.  Similar  definition  applies  to  Fj^(-Xj^,t). 
Equation  (40)  can  now  be  inverted  as 
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a^(2na)  +  x^,t)  =  -j 


1 

*  2 


1 

"  2 


‘a,(2»,t-i-t')dD,(x,..'*i) 


(42) 


Equation  (42)  can  be  written  in  a  compact  form  if  we  observe  that  a^(2n£D,t) 
and  Vj^(2naj,t)  vanish  for  t  <  0,  and  Dj^{±Xj^,t)  and  Fj^(±Xj^,t)  vanish  for 
t  <  ±  that  -hj^  <  Xj^  <  h^  by  Eq.  (18) .  We  have 

1  f' 


aj(2iK0.x,.t)  •■jJ  o^(2no,t-t')  d|Dj(-Xj,t')  .Dj(Xj,t')| 

'■■^I'^'^IO^ 

t 


1 

■"  2 


Vj^(2naj,t  -  t')  d|F^(-Xj^,t')  -  F^  (x^,t')| 


By  a  similar  argument,  we  obtain  from  Eq.  (38b), 

1  r' 


v^(2noj  +  x^,t)  =  jJ  Vj(2nco,t  -  t')  d|D^(-x^,t')  +  D^(x^.t')| 


C-h^/Ci^)- 


i  J  a^(2n(o,t-t')  d|E^(-x^.t')  -E^(x^,t')| 
(-hi/c^O^ 


(43a) 


(43b) 


When  material  1  is  elastic,  m^^  is  a  constant  and  ,  E^  and  F^^  are 


step  functions; 
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D^(±  x^.t) 


E^Ci  Xj,t) 


F^C±  x^.t) 


H  t  + 


'10 


.±tiL^lL\ 

"*1  V  ‘^10/ 


H  t  ?  - 


10 


Equations  (43)  then  reduce  to  (when  material  1  is  elastic) 

ajC2nc«.x^,t)  =  4  L(2na.,t  .  | 


(44) 


C45a) 


v^(2n4)>x^.t)  =  ^v^(2n«.t. 


(45b) 


This  is  nothing  more  than  the  familiar  characteristic  relation  for  one¬ 
dimensional  waves  in  a  homogeneous  elastic  medium.  As  such,  Eq.  (45)  may 
be  regarded  as  the  characteristic  relation  in  an  integral  form  for  one¬ 
dimensional  waves  in  a  homogeneous  viscoelastic  medium. 

•After  finding  the  stress  and  velocity  at  the  center  of  the  5th  layer 
in  Figs.  3  and  4,  we  use  Eq.  (45a)  to  find  the  stress  history  at  the  inter¬ 
face  between  the  4th  and  the  5th  layers  by  letting  n  =  4  and  x^^ 

The  result  is  shown  in  Fig.  6  along  with  the  exact  solution  by  the  ray  theory. 


7.  DISCUSSION  AND  CONCLUDING  REMARKS 

Several  analogies  between  the  solutions  of  transient  waves  in  a  finite 
layered  medium  and  in  a  finite  homogeneous  viscoelastic  medium  are  established. 
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The  analogies  between  the  solutions  apply  to  the  centers  of  the  alternate 
layers  in  the  layered  medium.  Solutions  at  points  other  than  the  centers 
of  the  alternate  layers  are  obtained  in  terms  of  the  solutions  at  the  centers 
of  the  alternate  layers  through  the  use  of  convolution  integrals  with  auxi¬ 
liary  functions  introduced  in  the  paper.  The  materials  in  the  individual 
layer  of  the  layered  medium  can  be  elastic  or  viscoelastic,  although  numeri¬ 
cal  examples  are  given  for  an  elastic  layered  medium. 

The  viscoelastic  analogies  derived  here  lead  us  to  an  exact  expression 
G(p),  which  is  the  Laplace  transform  of  the  relaxation  function  GCt)  for 
the  "equivalent"  homogeneous  viscoelastic  medium.  The  relaxation  function 
GCt)  is  obtained  by  numerically  inverting  its  Laplace  transform  using  the 
method  outlined  in  [13,14] .  The  method  used  in  [13,14]  tends  to  provide 
the  value  of  G(t)  only  at  a  finite  number  of  t’s  which  are  close  to  t  =  0. 

It  provides  poor  information  on  G(t)  for  large  t.  Fortunately,  our 
relaxation  function  approaches  rapidly  to  the  equilibrium  relaxation 
modulus  G^  so  that  the  inversion  obtained  by  the  method  in  [13,14]  is 
adequate  for  most  cases.  One  could  certainly  obtain  a  better  numerical 
inversion  of  a  Laplace  transform  by  using  other  techniques  such  as  the  fast 
Fourier  transform  [16]. 

Even  though  the  relaxation  function  G(t)  obtained  here  is  quite  crude, 
we  have  reached  an  excellent  agreement  between  the  solutions  by  the  visco¬ 
elastic  analogy  and  by  the  exact  ray  theory.  The  differences  in  the 
solutions  are  caused  mainly  from  the  numerical  Laplace  inversion  of  the 
relaxation  function  G(p)  and  the  auxiliary  functions  M(p)  and  Jj^(p).  For 
those  solutions  which  require  no  auxiliary  functions,  the  differences  in 
the  solutions  are  less  noticeable.  Some  of  the  auxiliary  functions,  namely, 
M,  R,  Lj^  and  defined  in  Eq.  (15)  can  be  determined  exactly  when  the 
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constituents  of  the  layered  medium  are  elastic.  With  the  use  of  e.xact 
auxiliary  functions,  the  differences  between  the  solutions  by  the  visco¬ 
elastic  analogy  and  by  the  ray  theory  can  be  made  smaller  [14], 

One  might  ask  the  relative  advantages  of  the  viscoelastic  analogy 
over  a  direct  numerical  computation  of  the  original  layered  problem. 

In  a  direct  computation  of  waves  in  the  original  layered  medium,  the 
calculations  are  feasible  when  the  boundary  conditions  are  constant  in 
tine  and  the  individual  layers  are  elastic.  Moreover,  keeping  track  of 
every  reflection  and  transmission  of  waves  at  the  interfaces  between  the 
layers  may  soon  exhaust  the  storage  capacity  of  the  computer,  not  to 
mention  the  computing  time.  The  situation  is  particularly  acute  when 
there  are  a  large  number  of  layers  involved  as  in  a  real  composite. 

These  shortcomings  are  not  present  in  the  theory  of  viscoelastic  analogy. 

One  possible  shortcoming  of  the  theory  of  viscoelastic  analogy  approach 
is  the  less  accurate  result  for  the  solution  at  points  other  than  the 
centers  of  the  layers.  This  shortcoming  is  not  serious  in  practical 
Implications  'vhen  the  individual  layers  are  very  thin  such  that  the 
solution  at  the  center  of  a  layer  and  at  other  points  in  the  layer  are 
practically  the  same. 

In  connection  with  the  present  work  on  transient  waves  in  a  finite 
layered  medium,  we  would  like  to  point  out  that  harmonic  waves  in  a  finite 
layered  medium  has  been  considered  by  Herrmann,  Beaupre  and  .■Vuld  [17].  In 
contrast  to  the  normal  stress  waves  studied  here,  the  harmonic  waves  con¬ 
sidered  in  [17]  are  horizontally  polarized  shear  waves.  However,  one  can 
see  from  the  analyses  presented  here  that  if  we  replace  the  normal  stress 
and  the  normal  displacement  by  the  shear  stress  and  the  transverse  displace¬ 
ment,  respectively,  the  analyses  presented  here  apply  to  the  transient 
shear  waves  in  finite  layered  media  as  well. 
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Viscoelastic  Relaxa 


Stress  at  the  center  of  the  Svh  layer  due  to  a  unit  step  stress  a])plied  at  the  center 
of  the  first  layer  while  the  center  of  the  14th  layer  is  fixed. 
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CHAPTER  7 

THE  EFFECTS  OF  DISPERSION  AND  DISSIPATION  ON  WAVE  PROPAGATION 
IN  VISCOELASTIC  LAYERED  COMPOSITES 

ABSTRACT 

In  Chapter  6,  stress  response  at  a  finite  distance  from  the  impact 
end  in  an  elastic  or  viscoelastic  layered  composite  is  studied.  In  this 
Chapter,  the  stress  response  at  a  large  distance  from  the  impact  end  of  a 
viscoelastic  composite  is  investigated.  If  the  distance  is  not  large  enough, 
the  stress  response  is  oscillatory  due  to  the  dispersive  nature  of  the  com¬ 
posite.  As  the  distance  increases,  the  dissipation  effect  of  the  visco¬ 
elastic  materials  becomes  pronounced  and  eventually  wipes  out  completely 
the  oscillatory  response.  The  transition  from  the  oscillatory  response  to 
the  monotonic  response  is  controlled  by  a  parameter  y  which  contains  (a) 
the  impedence  mismatch  of  the  composite  which  contributes  to  the  dispersion, 
(b)  the  dissipative  properties  of  the  viscoelastic  materials  and  (c)  the 
distance  traveled  by  the  wave. 

1.  INTRODUCTION 

Consider  a  semi-infinite  periodic  layered  composite  as  shown  in  Fig.  1 
in  which  each  period  2a)  consists  of  two  layers  of  homogeneous,  isotropic, 
linear  viscoelastic  materials.  The  thickness  of  individual  layers  are  2h^ 
(i=l,2)  where  the  subscripts  1  and  2  refer  to  materials  1  and  2,  respective¬ 
ly.  We  will  consider  plane  wave  propagation  in  the  direction  x  which 


is  normal  to  the  layers.  For  the  problem  considered  here,  the  surface  .x  =  0 
need  not  be  the  central  surface  of  the  first  layer.  We  will  assume,  however, 
that  the  first  layer  in  which  x  =  0  is  located  is  occupied  by  material  1. 

The  composite  is  initially  at  rest  and  at  time  t  =  0,  time-dependent, 
uniformly-distributed  normal  and  shear  stresses  are  applied  at  the  surface 
x=0.  Since  the  problem  considered  is  linear,  the  solutions  due  to  the 
applied  normal  stress  and  the  shear  stress  can  be  treated  separately.  The 
two  solutions  are  mathematically  identical.  Therefore,  we  will  consider 
only  the  solution  due  to  the  applied  normal  stress  at  x=0.  Moreover,  we 
will  assume  that  the  applied  normal  stress  at  x=0  is  the  Heaviside  unit 
step  function  in  time  t,  because  the  solution  for  a  more  general  applied 
normal  stress  can  be  obtained  by  a  linear  superposition. 

The  stress  response  at  a  position  x  which  is  sufficiently  large  can 
be  obtained  by  an  asymptotic  analysis.  When  both  layers  are  elastic,  the 
solution  can  be  expressed  in  terms  of  an  integral  of  an  Airy  function  [1,2]. 
The  stress,  as  a  function  of  time  t,  oscillates  around  the  Heaviside  step 
function.  When  one  or  both  layers  are  viscoelastic,  the  asymptotic  solution 
can  ce  expressed  in  terms  of  an  error  function  Cl  ,2-2.  The  stress  response  is 
no  longer  an  oscillatory  function  of  t,  but  a  monotonically  increasing 
function  -of  t  which  approaches  to  the  unit  stress  as  t  increases. 

Since  elastic  materials  are  special  cases  of  viscoelastic  materials, 
one  might  ask  how  a  monotonic  solution  becomes  an  oscillatory  solution  when 
the  viscoelastic  materials  become  elastic.  Alternately,  one  might  ask  what 
would  be  the  behavior  of  the  asymptotic  solution  if  the  relaxation  functions 
of  the  viscoelastic  materials  are  nearly  step  functions.  Clearly,  when  the 
position  X  is  not  large  enough,  the  dissipative  effect  of  the  viscoelastic 
materials  does  not  have  enough  time  to  prevail  and  the  stress  response  is 
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essentially  governed  by  the  dispersive  nature  of  the  composite  which  causes 
the  solution  to  be  oscillatory.  As  .x  increases,  the  dissipative  effect, 
no  matter  how  small,  becomes  prominent  and  dampens  the  dispersive  mechanism 
so  that  the  solution  is  non-oscillatory.  The  purpose  of  this  paper  is  to 
study  the  effects  of  the  dispersion,  dissipation  and  the  distance  of  wave 
propagation  have  on  the  wave  profile.  To  simplify  the  analysis,  we  will 
consider  only  solutions  at  x  =  2ojN’  where  .M  is  an  arbitrary  positive 
integer. 

It  should  be  pointed  out  that  a  similar  problem  was  studied  by  Sve  [1] 
for  two  special  viscoelastic  materials.  The  imaginary  part  of  the  wave 
number  for  the  special  materials  is  assumed  to  be  proportional  to  the  abso¬ 
lute  value  of  the  frequency  or  proportional  to  the  square  of  the  frequency. 
Hegemier  [3]  obtained  as.vmptotic  solutions  for  elastic  composites  as  well 
as  viscoelastic  composites.  However,  his  solutions  differ  from  that  obtained 
here  and  in  [2].  A  discussion  on  the  differences  will  be  given  later. 


2.  SOLUTION  FOR  x  =  2cjN 

The  equations  of  motion  and  the  continuity  of  the  displacement  are 
given  by 


3a. 

1 


(i  =  1.2)  (1) 


Ci  =  1.2)  (2) 


where  a^,  £^,  v^,  (i  =  l,2)  are  the  normal  stress,  normal  strain,  normal 
particle  velocity  and  mass  density,  respectively,  A  dot  stands  for  differen¬ 
tiation  with  respect  to  time  t.  The  initial  and  boundary  conditions  are 


a^(x,0')  =  v^(x,0')  =  E^(x,0')  =  0 


(i  =  1.2) 


(3) 
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a.(oo,t)  =  0 


(i  =  1,2: 


fij 


a^(O.t)  =  H(t] 


f5) 


where  H(t)  is  the  Heaviside  unit  step  function.  The  relation  between 
and  is  written  in  the  form  of  Stieltjes  integral 


cr.  (x.t)  =  g.  (t-f)  de^  (f) 


(i  =  1.2} 


(6) 


where  gj^(t}  are  the  relaxation  functions  of  the  viscoelastic  layers. 
Let  f(p)  be  the  Laplace  transform  of  f(t}; 


(P)  = 


f(t}e'P^  dt 


0 


(7) 


Equations  (1-6)  then  reduce  to 
d^c. 


—  =  k?  a. 


3x^  1  1 


(S) 


a.(=°,p}  =  0  , 


^1  "  V 


where 


(9) 


ki  =VpiP7g” 


(10) 


Since  is  periodic  in  x  with  periodicity  Zoi,  by  using  the  Floquet 

theory  [4]  the  solution  for  x  =  2a)N  where  .N  is  an  arbitrary  positive 
integer  is 

<^l(x,p)  =  -  e  ^11) 

where  ic  is  the  characteristic  exponent  given  by  (see  [2,5,6]) 

cosh  (2u)K)  =  9  cosh  (2kjh^  +  2k.,h2)  -  (9-1)  cosh  (2k^h^  -  2k.,hT)  (12) 
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9 


l=2'=i 

577 


Therefore,  the  solution  for  x  =  2(*)N  is 


ai(x.t)  =  ^ 


pt-Kx 


dp 


Br 


(13) 


(14) 


For  a  large  x,  the  main  contribution  to  the  Bromwich  contour  integral 
of  Eq.  (14)  appears  to  come  from  the  values  of  integrand  near  p  =  0. 

Hence  we  must  study  the  behavior  of  <  near  p  =  0  before  we  evaluate 
Eq.  (14)  for  large  x. 


3.  BEHAVIOR  OF  K  .NEAR  p  =  0 

For  most  viscoelastic  materials,  the  relaxation  function  g^(t), 

(i =  1,2)  is  a  moriotonically  decreasing  function  of  t.  Let  be  the 

value  of  g^(t)  at  t =  “.  For  most  viscoelastic  solids  g^^  is  non-zero. 
If  g^(p)  is  the  Laplace  transform  of  g^(t). 


Pg^Cp)  =  P 


gj^Ct)e'P^  dt 


[g- (t)  -  g- Je'P^  dt 


For  small  p,  e  =  1  -  pt  +  . . .  Hence 


Pg^  (P)  =  gi„(l  *  a.  p  -  a^T^p-^  +  . . . ) 


where 


a.  g .  = 

1*100 


[gi(t)  -  g.J  dt 


T.  = 


^  ^2ico 


[gi(t)  -  g.Jt  dt 


(15) 


(16) 


(17) 
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It  is  seen  that  a.g.  is  the  area  between  the  curve  g.  (tl  and  the 
horizontal  line  gj^(t)=g^ao  while  is  the  distance  of  the  centroid  of  this 

area  from  t =  0.  According  to  [7],  a^  provides  a  measure  of  the  "viscosity" 
of  the  viscoelastic  materials.  An  example  of  relaxation  function  which 
yields  Eq.  (16)  is  the  standard  linear  viscoelastic  solid 


(18) 


Using  Eqs.  (16),  (10)  and  (13),  the  right-hand  side  of  Eq.  (12)  can 
be  expanded  into  power  series  in  p.  If  we  assume  that,  for  small  p,  <  can 
be  e.xpressed  as 


0  2 

o<  =  P  -  JfP 


(19) 


and  use  of  this  to  expand  the  left-hand  side  of  Eq.  (12)  into  power  series  in 
p,  we  can  determine  the  constants  c^,  u  and  S  by  comparing  the  coefficients 
of  same  powers  of  p  on  both  sides  of  Eq.  (12).  After  a  lengthy  algebra, 
one  obtains 

<  =  gJP  .  (20) 


where 


(21) 


(22) 


(23) 
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We  see  that  p  and  are,  respectively,  the  effective  mass  density  and 

the  effective  equilibrium  modulus  of  the  composite. 

When  both  layers  are  elastic,  a^  =  0  and  hence  u  =  0.  Moreover,  only 
the  first  term  of  3  remains  and  3>0.  Notice  that  the  first  term  of  3 
is  proportional  to  the  difference  in  the  impedances  of  the  two  layers  and 
becomes  zero  when  the  difference  in  the  impedances  is  zero.  Since  the  dis¬ 
persive  nature  of  the  composite  comes  from  the  impedance  mismatch,  the  first 
term  of  3  is  responsible  for  the  oscillatory  nature  of  the  stress  response. 

When  one  or  both  of  the  layers  are  viscoelastic,  u  is  positive  and 
non- zero  while  p  can  be  positive,  negative  or  zero.  Not  only  is  u 
responsible  for  the  dissipative  nature  of  the  stress  response,  the  second 
part  of  3  is  also  responsible  for  the  dissipation. 

The  case  when  both  u  and  3  vanish  will  not  be  considered  here. 

4.  ASYMPTOTIC  SOLUTIONS 

From  Eqs.  (19)  and  (14),  we  have 

Br 

We  will  assume  that  .x  is  sufficiently  large  that  the  terms  denoted  by  the 
dots  can  be  ignored.  We  will  also  assume  that  3^0.  The  case  3  =  0  will 
be  discussed  later.  Let 


1/3 


(25) 


217 


Equation  (24)  then  takes  the  form 


o'iy.T) 


1 

2iTi 


i  +  ±5  p' 


dp 


(26) 


Br 


where  the  subscript  1  of  a  has  been  omitted  and  the  +  sign  is  for  B  >  0 
and  -  sign  for  B<0. 

By  taking  the  Bromwich  contour  as  shown  in  Fig.  2,  it  is  not 

difficult  to  show  that 

a'"(Y,T)  +  a'(Y,-x)  =  1  (27) 

We  will  therefore  consider  only  the  case  S>0  and  hence  the  integral 

—  e 
P 

Br 

Using  the  identity 


r  I  Tp  +  YP'^-^iP^ 
p  ®  ^  dp 


(28) 


i  e^P  =  i  >  f  e^P 


e  ds 


(29) 


the  integral  in  Eq.  (28)  can  be  divided  into  two  parts: 


a(Y,T)  =  +  I.,  , 


(50) 


I,  = 


1 

2TTi 


1  YP^-^ip^ 

-  e  dp 


Br 


(31) 


4  I 

ds  ! 

1 

i  1 

27:i 

E 

2  1  1 
C  sp  +  YP'^+JP* 

'  e  dp 

Br 


(32) 


Notice  that  Ij  ■  c(y,0)  and  hence  is  the  magnitude  of  stress  at  t =» x/c^.  By 
taking  the  Bromwich  contour  L,  of  Fig.  2,  one  obtains 
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I,  = 


i  +  i 

3  *  t: 


i  .  ^ 

3  *  27r 


ie  2 

r 


1  2^3 


yr 


sin(-^  yr  j  dr 


(33) 


3-1/3  r(|),.3‘/^r(l)  If.,.,  j 


where  r(x)  is  the  Ganuna  function.  If  y  is  very  large,  we  take  the 
Bromwich  contour  and  obtain 


0 

1  6 


(34) 


where 


2  12/F  I  144 

/2c  B^1/2 


l-#r6^>... 


5  =  Y 


-3/ 


XU 


(35) 


We  now  turn  to  the  integral  l2-  5/  replacing  the  variable  p  by 


p  =  z  -  y  , 


(36) 


Equation  (32)  can  be  written  as 


2  3 

-sy  +  ^y 

=  I  e  3  ^  ds 
0 


J_f 


2ni 


I 

Br 


dz 


(37) 


or 


^2 


2  3 

-sY+-y  - 

e  Ai  (-S  +  y  )  ds 


(38) 


where  the  Airy  function  is  defined  as  [8] 


Ai(s)  = 


1 

2TTi 


1  3 

-SZ  +  -z-Z 

e  dz 


(39) 
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IT 


COS 


(sr.ir’) 


dr 


(39) 

(Cont 'd) 


Tw  extreme  cases  of  y  =  0  and  y  =  ®  have  been  studied  in  the 
literature.  Before  we  evaluate  a  for  arbitrary  y,  we  will  obtain  these 
two  extreme  cases  from  Eq.  (30). 


(1)  y  =  0 

For  elastic  composites,  u=0  and  hence  y  =  0.  Equations  (30),  (33) 
and  (58)  then  yield 


a(0,T) 


3 


T 

Ai(-s)  ds 

0 


(40) 


This  is  precisely  the  asymptotic  solution  obtained  in  [1,2].  The  stress  a 
is  an  oscillatory  function  of  t  (see  Fig.  3). 


(2)  y  =  «> 

For  viscoelastic  composites,  u 0  and  S  of  Eq.  (24)  may  or  may 
not  be  zero.  In  [2]  the  term  containing  S  was  ignored.  This  is  equiva¬ 
lent  to  assuming  that  3  =  0  and  hence  '(■=«>.  For  a  very  large  y,  the 
•Airy  function  has  the  e.xpression  [8]; 

-  1  1  -  4- 

Use  of  this  expression  in  Eq.  (38)  results  in 


0 


1 


where 


Therefore,  when  0  =  0,  (i.e.,  y  =  “  or  6=0),  Eqs.  (34),  (42)  and  (30)  yield 

the  following  asymptotic  solution  obtained  in  [1,2]  ; 

a  =  3  ^  1  +  erf(T  /2)J-  (44) 

The  stress  a  is  a  monoton i ca 1 1 y  increasing  function  of  ,  Fig.  4. 


5.  NUMERICAL  RESULTS  AND  DISCUSSION 

For  an  arbitrary  y.  The  stress  o  as  a  function  of  t  may  be 
obtained  from  Eq.  (30)  where  I^  is  given  by  Eq.  (33)  or  (34)  and  I^  is 
given  by  Eq.  (38).  Since  both  I^  and  require  a  numerical  integra¬ 

tion,  it  might  be  simpler  to  evaluate  a  directly  from  Eq.  (28).  If  we 
take  L^  of  Fig.  2  as  the  Bromwich  contour,  Eq.  (28)  reduces  to 


o(y,T)  = 


2 

1  y(Tr-yr^ 
—  e 
r 


Sin|^-y 


(xr  +yr^) 


For  the  contour  L^,  we  have 


o(y,T)  'r 

0 

where  x*  and  5  are  defined  in  Eqs.  (43)  and  (35).  Notice  that 

X  =  X*  when  y  =  6  =  1  .  (47) 

.J.2 

Notice  also  that  because  of  the  factor  (l/r)e  in  the  integrand  of 

Eq.  (46),  the  absolute  value  of  the  integrand  diminishes  rapidly  as  r 

-r*  -  -3 

increases.  For  instance,  at  r  =  2,  (l/r)e  =9x10  and  at  r=3. 
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/  - -  * 5 

(l/r)e  =4x10  .  Therefore,  the  infinite  integral  can  be  replaced  by  an 

integral  of  finite  interval  say  0  <  r  <  3.  A  similar  argument  applies 
to  the  integral  in  Eq.  (45). 

Equation  (45)  is  used  to  calculate  a  for  Y  =  0,  0.1,  0.3,  0.6  and 
1.0.  The  results  are  shown  in  Fig.  3.  Equation  (46)  is  used  to  calculate 
a  for  Y=1  and  Y  =  (i.e.,  6  =  1  and  6  =  0),  Fig.  4.  We  see  that  the 

stress  response  differs  very  little  for  Y=1  and  y  =  “• 

The  e.xample  of  stress  response  at  the  30th  layer  considered  in  [2] 
has  a  negative  value  of  Q  and  y=0.58.  On  the  other  hand,  the  example 
considered  on  p.  110  of  [3]  has  a  positive  value  of  8  and  '{  =  0.68. 

For  a  given  viscoelastic  composite,  u  and  8  are  known  and  fixed. 
Y  then  depends  on  x  and  increases  as  x  increases.  We  see  from  Fig.  3 
that  the  oscillatory  nature  of  the  stress  diminishes  as  y  increases. 
Since  for  y8  1  the  oscillation  is  practically  non-existence,  we  may  say 
that  for 


the  stress  response  is  monotonic. 

The  asymptotic  solution  for  viscoelastic  composite  derived  in  [3]  is 
different  from  Eq.  (24).  Using  the  notations  of  Eqs.  (20-22),  the  asympto¬ 
tic  solution  derived  in  [3]  is  based  on  the  equation 


a(x,t) 


1 

2iTi 

Br 


1 

—  exp 
P 


tp 


2c 


dp 


(49) 


If  we  expand  the  last  term  in  the  exponent  into  a  power  series  in  p  and 
ignore  the  terms  of  order  higher  than  p^,  Eq.  (49)  is  identical  to 
Eq.  (24).  We  are  able  to  verify  that  J  in  [3]  is  identical  to  the  one 
obtained  in  Eq.  (21).  However,  8  in  [3]  appears  to  be  different  from  the 
expression  in  Eq.  (22). 
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6.  CONCLUSION 

The  parameter  y  defined  in  Eq.  (25)  consists  of  the  variables  u,  S 
and  X.  The  dissipative  nature  of  the  viscoelastic  material  is  represented 
by  u  and  a  part  of  6,  while  the  remaining  part  of  8  represents  the 
dispersive  nature  of  the  composite.  The  distance  traveled  by  the  wave  is 
represented  by  x.  Thus  y  contains  the  influences  on  the  wave  profile 
due  to  dissipation,  dispersion  and  the  distance  traveled  by  the  wave.  With 
y  determined  from  Eq.  (25),  Figs.  3  and  4  provide  the  wa.ve  response. 
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CHAPTER  8 
CONCLUSIONS 


Static  and  dynamic  analyses  of  multilayer  composite  laminates  and/or 
pilot  studies  leading  toward  such  analyses  have  been  presented  in  this  report. 
Detailed  summaries  and  conclusions  have  been  included  with  each  chapter.  The 
present  chapter  is  intended  to  briefly  summarize  these  conclusions  and,  where 
possible,  to  make  suggestions  for  further  investigations. 

Three  approaches  were  investigated  in  Chapter  2  for  the  finite-element 
analysis  of  edge-effects  in  symmetric  laminates  under  prescribed  inplane  strain. 
Although  these  approaches  produce  essentially  identical  results  away  from  the 
free  edge,  some  differences  are  noted  in  the  magnitude,  and  in  selected  cases 
the  form,  of  the  stress  distributions  near  the  free-edge.  Until  further  analytic 
results  are  available  it  is  not  possible  to  conclude  which  of  the  strain  contin¬ 
uity  or  traction-free-edge  approaches  is  the  more  accurate.  Some  insight  into 
this  comparison  may  be  obtained  by  combination  of  the  stress  predictions  with 
appropriate  failure  criteria;  the  better  analysis  should  be  capable  of  consis¬ 
tent  predictions  of  laminate  first-failure  in  comparison  with  experimental 
results. 

Of  the  8-node  single  layer  plate  bending  elements  with  a  straight  traction- 
free  edge,  developed  in  Chapter  3,  the  best  element  is  identified  as  one  based 
on  a  21^  stress-field.  This  element  produced  displacement  and  stress  predic¬ 
tions  which  were,  in  general,  superior  to  all  other  elements  tested.  Although 
the  superiority  of  this  element  over  a  standard  hybrid-stress  plate  element 
is  not  evident  in  all  cases,  it  is  expected  that  the  superiority  will  be 
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apparent  in  examples  where  free-edge  effects  are  dominant,  such  as  in  multi¬ 
layer  laminates.  In  the  development  of  an  8-node  multilayer  plate  element, 
a  layer  stress  field  analagous  to  that  used  in  the  21 B  element  is  recommend¬ 
ed. 

In  Chapter  4,  alternate  initial-stress  approaches  based  on  the  hybrid- 
stress  model  were  compared  for  the  elastic-plastic  analysis  of  single  layer 
plates.  Results  obtained  for  selected  example  problems  suggest  that  hybrid 
functional  II,  coupled  with  iteration  scheme  A,  is  the  preferred  approach. 
This  approach  should  therefore  be  used  in  extension  to  include  material  non¬ 
linear  effects  in  multilayer  laminated  plates. 

Chapter  5  provides  a  means  for  analysing  the  natture  of  singularities 
at  the  free  edge  of  a  composite  whose  individual  layers  are  anisotropic. 

The  results  obtained  can  be  used  to  formulate  a  hybrid-stress  singularity 
element  for  the  free-edge  point. 

The  two  dynamic  problems  studied  in  Chapters  6  and  7  contributed  signif¬ 
icantly  to  the  literature  of  dynamic  response  of  composites.  A  theory  of 
viscoelastic  analogy  presented  in  Chapter  6  offers  a  reliable  way  to  pre¬ 
dict  the  transient  response  of  a  layered  composite  at  finite  times.  More¬ 
over  the  dimension  of  the  layered  composite  is  finite.  This  problem  was 
not  solved  successfully  before. 

The  other  dynamic  problem,  namely  the  transient  response  of  a  visco¬ 
elastic  layered  composite  studied  in  Chapter  7,  sheds  much  light  on  the 
nature  of  interaction  between  the  effects  of  dispersion,  dissipation  and 
the  distance  traveled  by  the  wave.  The  results  show  how  an  oscillatory 
wave  approaches  a  monotonic  wave  as  the  distance  traveled  by  the  wave 
increases. 
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